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Boundary-value  Problems  II:  C^-cubics  with  Detailed  Results 

by 

Andrew  J.  Martin  and  James  W.  Daniel 

Abstract 

The  methodology  is  very  briefly  described  and  then  numerical  results 
are  presented  for  an  implementation  of  a smooth-cubic-spline-collocation  pro- 
cedure accelerated  via  iterated  deferred  corrections  to  obtain  approximate 
solutions,  accurate  to  a prescribed  tolerance,  of  two-point  boundary-value  prob- 
lems for  second-order  scalar  ordinary  differential  equations.  The  results  are 
similar  to  those  obtained  via  a less  generally  applicable  finite-difference- 
oriented  code  described  elsewhere  which  competes  well  with  the  best  codes  avail- 
able . 

1 . Introduc  tion 

The  primary  methods  for  the  approximate  solution  of  nonlinear  two-point 
boundary-value  problems  can  be  roughly  classified  into  four  types:  (1)  Rayleigh- 

Ritz-Galerkin  and  collocation  methods  [deBoor-Swartz  (1973),  Douglas-Dupont  (1974). 
Russell  (1974),  Russell- Shampine  (1972),  ut  cetera],  (2)  shooting  methods  [Bulirsch 
et  alia  (1975),  Keller  (1968),  Scott-Watts  (1975),  et  ce tera ] , and  (4)  standard 
finite-difference  methods  [Keller  (1968,  1974,  1975),  Lentini-Pereyra  (1974. 

1975a,  1975b),  Pereyra  (1968.1973).  et  cetera ] . Despite  the  dominance  of  the 
Rayleigh-Ritz-Galerkin  and  collocation  methods  in  the  more  theoretical  literature, 
there  seems  to  be  a consensus  [Aziz  (1975) 1 that  the  most  efficient  practical 
procedures  are  the  shooting,  imbedding,  and  difference  methods,  although  to  be 


competitive  the  difference  methods  must  be  implemented  with  an  acceleration 
technique  such  as  extrapolation  [Joyce  (1971),  Keller  (1974,  1975)]  or  deferred 
correction  [Fox  (1947,  1962),  Joyce  (1971),  Lentini-Pereyra  (1974,  1975a,  1975b), 
Pereyra  (1968,  1973),  Daniel-Martin  (1975,  1977)].  To  make  the  collocation  methods 
also  competitive  it  has  been  proposed  [Russell  (1974),  Daniel  (1974)]  that  they 
be  implemented  with  acceleration  techniques  much  as  are  difference  methods. 

In  [Daniel  (1974)],  the  second  of  the  present  authors  described  how  de- 
ferred correction  could  be  used  with  various  spline-collocation  methods,  although 
no  results  were  given  for  any  actual  implementation.  The  present  paper  summarizes 
the  numerical  results  obtained  with  a careful  implementation  of  deferred  correction 
applied  to  one  of  the  spline-collocation  methods  proposed  in  that  paper,  namely 
the  so-called  extrapolated  collocation  method  [Daniel- Swartz  (1973),  Hill  (1973)1 
which  uses  twice-continuously-differentiable  cubic  spline  functions  and  was  moti- 
vated by  the  work  of  [Fyfe  (1969)].  The  code  SPLIDC  (spline  _iterated  deferred 
correction)  implementing  this  is  a modification  of  the  code  NUMIDC  (Numerov ' s 
method  with  iterated  deferred  correction)  developed  by  the  authors  and  shown  in 
[Daniel-Martin  ( 19 75 > 19  77)]  to  be  quite  competitive  with  the  best  available  codes 
on  appropriate  classes  of  problems.  The  present  paper  compares  the  behavior  of 
SPLIDC  with  that  of  NUMIDC  on  the  test  problems  earlier  used  for  NUMIDC,  namely 
a set  of  some  seventeen  second-order  problems  not  involving  the  first  derivative 


of  the  unknown  function  (so  that  Numerov's  method  will  be  fourth-order  accurate): 
SPLIDC  often  follows  the  same  computational  paths  as  did  NUMIDC  and  is  about  as 
reliable,  although  SPLIDC  involves  more  effort  in  terms  of  function  evaluations 
and  much  more  total  time  because  of  overhead  in  solving  systems  of  linear  equations. 
We  also  present  the  results  of  using  SPLIDC  to  solve  some  problems  explicitly 
involving  the  first  derivative  in  the  differential  equation,  problems  for  which 
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Nuraerov's  method  drops  to  second-order  accuracy  (so  NIJMIDC  is  not  applicable) 
while  extrapolated  collocation  remains  at  fourth-order  accuracy  (so  SPLIDC  is 
applicable)  . Although  SPLIDC  was  designed  to  handle  the  involvement  of  the  first 
derivative  in  the  boundary  conditions,  we  have  observed  computationally  some 
puzzling  behavior  and  significantly  greater  numerical  difficulties  in  this  case; 
we  therefore  do  not  allow  such  boundary  conditions  in  SPLIDC  and  we  will  devote 
considerable  effort  to  resolving  these  difficulties  as  we  proceed  to  study  the 
use  of  different  spline  spaces  (see  the  following  paragraph) . 

Overall,  our  results  serve  as  a feasibility  study  for  using  deferred  cor- 
rection with  spline-collocation  methods;  they  demonstrate  that  the  methods  proposed 
in  [Daniel  (1974)]  can  in  fact  work  in  practice.  What  is  needed  next  is  to  de- 
termine what  spaces  of  spline  functions  can  be  used  most  efficiently  as  the  basis 
of  such  methods.  We  are  pursuing  such  questions  now  and  will  report  later  on  the 
results . 

In  the  next  section  of  this  paper  we  describe  briefly  the  method  imple- 
mented by  SPLIDC;  for  more  detail  the  reader  is  referred  to  [Daniel  (1974)].  In 
Section  3 we  briefly  describe  how  SPLIDC  implements  the  method  of  Section  2;  since 
the  structure  of  SPLIDC  is  a slight  modification  of  that  of  NUMIDC,  the  reader  is 
referred  to  [Daniel-Martin  (1975,  1977)]  and  to  the  code  for  SPLIDC  in  Section  8 


2 . The  Me  tjbod  to  Be  Imp  lemon  ted 

A detailed  description  of  iterated  deferred  correction  used  with  extra- 
polated collocation  may  be  found  in  Section  4 of  [Daniel  (1974) ];  specifically, 
the  material  from  the  paragraph  before  Equation  4.6  through  the  paragraph  follow- 
ing Equation  4.10  explains  the  process  used.  Here  we  will  be  very  brief.  We 
consider  the  differential  equation 

(2.1)  y"(t)  f(t,y(t),y'(t))  for  a < t < b 

subject  to  boundary  conditions 

(2.2)  y(a)  yQ  given,  y(b)  yp  given, 

and  we  consider  approximating  y by  a twice-continuously-dif ferentiable  cubic 
spline  Z on  a uniform  mesh 

(2.1)  a t.  < t,  < . . . < t b with  t.  , - t.  h 7;  for  0 < i < N-l; 

01  N i ‘ 1 i N — — 


thus  Z is  given  by  a cubic  polynomial  on  each  piece  (t^,t^  j)  for  0 < i < N- 1 . 

We  let  S be  another  cubic  spline  which  interpolates  the  solution  y to  Equation  2.1 
and  2.2  at  the  mesh  points  of  Equation  2.3  and  which  satisfies  some  additional 
special  boundary  conditions.  It  follows  [Daniel  (1974)]  that  S satisfies  differ- 
ential equations  similar  to  that  in  Equation  2.1;  for  example,  writing  g^  for 

g(t.),  we  have  (for  certain  parameters  e.  and  d .) 
i J J 

(2.4)  S"  t -^(SV  . -2S"  ■ S'.'  ) - q yf‘’j  ‘-‘e.  f(t  ,S  ,S!  y +5)  h2J  ‘4d  .) 

l 12'  i-l  i i'l  j 0 1 J 111  j o 1 J 

-i  0(h‘5q'1) 


for  i c i < N-l,  with  similar  equations  for  i 0 and  i 


N. 


A sequence  of 


approximating  splines  Z is  obtained  essentially  by  using  successively  more  terms 


(higher  values  of  q)  in  liquation  2.4,  yielding  approximations  accurate  at  the 

4 8 12 

mesh  points  to  successively  higher  orders:  9(h  > . 9(h  ),  0(h  ~) , et  cetera. 


More  specifically,  the  zero-th  leve  1 spline  , with  values  Z^  ^ at 


the  mesh  points,  satisfies  Z 


0,0 


y0  and  Z0.N 


y , and  solves 
r 


Zi'.l  1 Ts(2i',i-l-aZ0.liZ0.:L.l)  £tti’Z0.1-Zd.i1 

4 

witli  similar  equations  at  i 0 and  i N;  it  turns  out  that  Zq  S • 6(h  ) so 
4 4 

that  Zq  y i 0(h)  and  Z^  . y^  t 0(h  ) . The  deferred  correction  step  leads 

from  Z„  to  Z,,  the  solution  at  level  one,  via  requiring  Z,  „ y„ , Z,  y„,  and 

U 1 1,0  0 1 , N r 

<*•»  zu  ■ e<vzi,i'z;.i 


j 0 


for  1 < i < N-l,  with  similar  equations  for  i 0 and  i N;  denotes  an 

h ^ ^-accurate  difference  approximation  to  y^-^  computed  using  the  values  Z^  ^ 

8 ' u 

rather  than  the  values  y^.  It  turns  out  that  Z^  S + 0(h  ) so  that  Z^  y < CHI'  ) 


and  Z^  ^ yi  )•  Similarly,  further  deferred  correction  steps  lead  to  sue- 

4m  • 4 

cessively  higher  level  solutions  Z,,Z. satisfying  Z S + ©(h  ),  so  that 

tl  J m 

4 4m  * 4 

y * 9(h  ),  Z ^ y^  ) • It  is  also  possible  easily  to  compute 

asymptotically  correct  estimates  of  the  level-m  errors  E . y.  - Z . essen- 

m, l i m, i 

tially  by  performing  one  step  of  Newton's  method  for  finding  Z , . - Z 

m+1, i m, i 

The  simplest  implementation  of  iterated  deferred  correction  for  extrapolated 
collocation,  given  a user-provided  error  tolerance  TOL,  would  thus  somehow  select 
an  initial  h ^ and  start  computing  Z^.E^.Z  . . . on  this  mesh,  ideally  until 


'I 
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Z ( via  K ) meets  the  error  criterion;  it"  this  does  not  succeed  or  at  least  pro- 
m ra 

gross  satisfactorily,  the  code  could  increase  N (for  example,  by  doubling)  and 
start  over  with  the  resulting  new  grid.  This  is  essentially  how  Numerov's  (or 
Cowell's)  method  was  implemented  in  [Pereyra  (1973)].  Since  our  code  NUMIDC 


[Daniel-Martin  (1975,  1977)]  for  that  method  incorporated  many  improvements  in 
overall  algorithm  design,  some  of  which  were  suggested  in  [Lentini-Pereyra  (1974, 
1975a,  1975b) |,  we  likewise  incorporated  these  improvements  for  extrapolated 
collocation  in  our  code  SPLIDC.  We  now  turn  to  that  code. 


3.  The  Implementation  of  the  Method 

A detailed  description  of  the  practical  modifications  in  the  basic  ap- 
proach to  iterated  deferred  correction  for  Numerov's  method  may  be  found  in  Sec- 
tions 3 turough  6 of  [Daniel-Martin  (1977)]  and,  in  somewhat  more  detail,  in 
Sections  3 through  7 of  [Daniel-Martin  (1975)].  Our  code  SPLIDC  implementing 
iterated  deferred  correction  for  extrapolated  collocation  follows  almost  exactly 
that  same  structure  as  outlined  for  NUMIDC  in  Figure  7.1  of  [Daniel-Martin  (1977)]. 
Some  slight  modifications  in  structure  were  required  in  going  from  NUMIDC  to 
SPLIDC  to  account  for  the  differences  in  the  two  basic  methods,  the  fact  that 
SPLIDC  allows  y'  to  appear  in  the  differential  equation,  et  cetera.  We  still 
allow  at  most  only  256  intervals  in  the  mesh.  A thoroughly  commented  code  for 
SPLIDC  may  bo  found  in  Section  8 and  is,  of  course,  the  ultimate  authority  on 
the  structure  of  our  algorithm's  implementation.  The  fundamental  improvements  in 
SPLIDC,  as  in  NUMIDC,  and  as  compared  with  the  simplest  implementation  sketched 
at  the  end  of  the  preceding  section,  are:  (1)  avoiding  useless  deferred-correction 

steps  at  high  levels  by  checking  to  see  that  the  numerical  solutions  are  behaving 
at  the  proper  order  asymptotically  in  h before  increasing  the  level;  (2)  avoiding 


level  on  the  new  mesh;  (1)  solving  the  nonlinear  equations  more  efficiently  and 
reliably;  and  (4)  improving  the  control  of  both  discretization  errors  and  round- 
ing errors. 

4 . Summary  of  Numerical  Results 

In  this  section  we  describe  the  overall  results  of  our  experience  with 
SPLIDC  and  give  some  examples  to  illustrate  specific  features  of  the  implementa- 
tion. Data  on  the  final  accuracy  and  the  cost  to  achieve  it  for  each  individual 
problem  are  presented  in  Section  7,  so  we  merely  summarize  here;  detailed  infor- 
mation on  the  paths  followed  during  the  computation,  the  number  of  Newton  itera- 
tions, et  cetera,  is  available  from  the  authors.  All  problems  were  solved  by 
our  code  SPLIDC  on  the  CDC  6600  at  The  1'niversity  of  Texas  at  Austin  using  the  RUN 
compiler;  the  execution  times  given  exclude  output  but  are  still  somewhat  inde- 
terminate (to.  say.  within  ’-10%)  because  timing  involves  a system  request.  In 
each  case  our  codes  attempt  to  approximate  the  true  solution  y by  a spline  7. 

within  a given  tolerance  TOL  in  the  sense  that  we  desire  max  |y.  -z-  I - TOL*  max  |Z. 

O-'i'-N  1 1 0-i-N  1 

Our  first  table.  Table  4.1.  summarizes  the  performance  of  SPLIDC  on 
thirteen  non-pathological  problems  from  [Daniel-Martin  (1975)],  all  but  three  of 
which  involve  nonlinearities  in  y (polynomial,  exponential,  or  logarithmic) , but 
none  of  which  explicitly  involve  y’ . SPLIDC  is  provided  with  a logical  switch  so 
that  the  user  can  inform  it  if  y1  is  absent  from  the  equation,  allowing  consider- 
able computational  savings  as  should  be  obvious  from  Equations  2.4  and  2.5.  We 
present  for  each  ot  four  tolerances  the  total  time,  total  number  of  evaluations 
of  f,  and  total  number  of  evaluations  of  x—  , needed  to  solve  all  thirteen  ot  the 


s 


problems;  we  also  list,  for  each  tolerance,  the  total  number  of  occurrences  of 
UNATTAINABLE  ANNOUNCED  (when  the  code  announces  that  it  cannot  attain  the  re- 
quested tolerance)  and  of  FAILURE  (when  the  code  incorrectly  announces  that  it 
has  attained  the  requested  tolerance).  This  summary  appears  in  Table  4.1;  for 
comparison,  we  present  the  analogous  data  [Daniel-Martin  (1977)]  for  the  finite- 
difference  code  NUMIDC  in  Table  4.2. 


Table  4 , 1 Performance  Summary  for  SPLIDC  on  Thirteen  Problems 


Table  4.2  Performance  Summary  for  NUMIDC  on  Thirteen  Problems 


TOLERANCE 

TIME (SECS.) 

f EVALUATIONS 

d£ 

EVALUATIONS 

UNATTAINABLE 

ANNOUNCED 

FAILURE 

io' 3 

0. 39b 

976 

946 

0 

0 

to'6 

1.156 

2674 

18  66 

0 

0 

io'9 

2.  341 

5264 

2762 

0 

0 

io-ly 

4 . 94  7 

10  623 

5244 

- - 

2 

0 

I;  mn  Table  4.1  we 
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can  see  that  iterated  deferred  correction  for  our  spline- 
collocation  method  as  implemented  in  SPLIDC  does  in  fact  work,  that  is,  it  does 
solve  most  of  the  test  problems;  this  provides  the  practical  justification  for 
[Daniel  (1974j 1 . Hy  comparing  'Table  4.1  with  Table  4.2,  we  also  see,  however, 
that  in  several  respects  SPLIDC  does  not  perform  as  well  on  these  problems  as 
does  NUMIDC: 

(1)  The  two  FAILURES  recorded  for  SPLIDC  occurred  because  the  actual  error 

was  four  to  five  times  the  estimated  error,  which  in  turn  barely  met  the  toler- 

-12  -12 

ance;  thus  an  accuracy  of  only  about  3 x 10  rather  than  10  was  achieved.  It 
should  be  noted  that  NUMIDC  similarly  underestimates  the  error  on  these  thirteen 
problems  about  10%  of  the  time,  but  by  coincidence  this  never  caused  FAILURE  in 
the  test  cases. 

(2)  The  much  greater  time  required  by  SPLIDC  than  by  NUMIDC  is  only 
slightly  due  to  the  greater  number  of  function  evaluations;  it  is  primarily  caused 
by  the  greater  overhead  involved  in  solving  the  more  complex  systems  of  equations 
involved  in  SPLIDC  than  in  NUMIDC. 

(3)  The  larger  number  of  function  evaluations  is  partially  due  to  the 
desire  to  give  this  program  the  flexibility  to  handle  y'  in  the  differential  equa- 
tion or  the  boundary  conditions.  The  current  version  reevaluates  the  function 
and  derivative  at  the  endpoints;  while  this  is  necessary'  if  y1  is  present  in  the 
function  or  the  boundary  conditions,  it  is  not  necessary  for  the  problems  used  to 
generate  Tables  4.1  and  4.2  and  could  rather  easily  be  eliminated.  If  these 
evaluations  were  eliminated  the  corresponding  entries  for  SPLIDC  would  be  as 


shown  in  Table  4.1a. 


If  we  make  this  modification,  then  for  the  47  cases  in  which  both  NUMIDC  and 
SPLIDC  reach  the  tolerance,  SPLIDC  requires  no  more  evaluations  than  does  NUMIDC 
in  11  cases,  or  6670  of  the  time. 

Thus,  it  appears  that  if  we  want  a deferred  correction  version  of  a spline 
collocation  process  to  compete  with  NUMIDC  for  problems  in  which  y'  is  not  ex- 
plicitly present,  we  must  use  a basic  collocation  process  that  is  more  efficient 
than  is  extrapolated  collocation.  We  are  presently  examining  the  many  spline- 
collocation  methods  available  with  the  goal  of  finding  such  efficient  processes. 

Recall  now  that  SPLIDC,  as  opposed  to  NUMIDC,  is  able  to  handle  equations 
explicitly  depending  on  y'.  Table  4.3  summarizes  the  behavior  of  SPLIDC  on  a 
set  of  six  such  test  problems  in  Section  7. 


TABLE  4,3  Performance  Summary  for  SPLIDC  on  Six  Problems  Involving  y' 


Tolerance 

Time  (secs.) 

f Eva luations 

df  df 

T-  or  "TT 

oy  2y 

Evaluations 

Unattainable 

Announced 

Fai lure 

10  ' 5 

0.67! 

693 

62  7 

0 

1 

10' 6 

2.0  75 

2431 

1711 

0 

1 

io'9 

i.954 

5219 

2421 

1 

0 

io  ’ l‘J 

5.964 

8957 

3240 

* 

2 

U-_ 

Again  we  see  chat  iterated  deferred  correction  of  extrapolated  collocation 

as  implemented  in  SPLIDC  does  work  to  solve  problems  involving  y'  explicitly  in 

the  equation,  thus  providing  the  computational  justification  lacking  in  [Daniel 

(1974)].  The  one  FAILURE  at  TOL  10  ( is  relatively  minor  in  that  a tolerance 

of  2x10  was  actually  achieved;  the  FAILURE  was  caused  by  the  fact  that  the 

true  error  was  about  two  to  three  times  the  estimated  error.  The  FAILURE  at 

TOL  10  was  negligible  since  the  tolerance  actually  obtained  was  1.0018x10 

again  caused  by  the  fact  that  the  true  error  was  larger  (by  a factor  of  less  than 

- 12 

1.33)  than  its  estimate.  One  of  the  FAILURES  at  TOL  = 10  was  slightly  sig- 

-12 

nificant  in  that  a tolerance  of  only  6.2x10  was  achieved,  due  to  true  error 

being  about  35  times  estimated  error;  the  other  FAILURE  was  minor,  since 
- 12 

1.3x10  was  actually  achieved. 

In  order  to  demonstrate  that  the  deferred  correction  process  is  actually 
working  and  that  SPLIDC  does  not  simply  remain  at  level  zero  with  fourth-oryder 
accuracy  from  the  basic  method,  we  present  now  some  examples  of  the  computational 
paths  followed  by  SPLIDC. 

-3 

At  a large  tolerance  like  10  , as  in  NUMIDC,  usually  the  tolerance  is  met 
at  level  zero  after  changing  the  mesh  spacing  one  or  two  times.  On  some  more 
difficult  problems,  correction  steps  may  be  performed.  This  is  exemplified  for 

-3 

TOL  10  in  Problem  #6  of  [Daniel-Martin  (1975)]: 


y y * y -i  e 


3 sin  4t rx  I 2 2.  . . 

-*  e [lorr  (cos  47TX  - sin  47tx)  - e 


2s  in  4ttx 


- lj  for  0 < x <■  1, 


y(0)  y(l)  1 with  solution  y(x)  e 


sin  4ttx 


The  performance  of  SPLIDC  is  described  in  Table  4.4. 
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TABU-:  4.4  PROBU;M  ?r6  WITH  TOL  10 ” 1 , SO  THAT  THE  MAXIMUM  ABSOLUTE  ERROR  MIST 

NOT  EXCEED  2.  7x  10  " 3 


An  entry  "a  t b(c)"  means  that  estimated  absolute  discretization  error  is  a,  a 

bound  on  the  error  in  solving  the  nonlinear  equations  G(z)  = 0 is  b,  while 

c is  the  actual  total  error.  It  is  a+b  that  estimates  the  total  error  and  is 

-3 

tested  against  2.7x10  = YNORMxTOL.  "NEWT  = p means  that  iterates  Zq,...,Z 

were  generated  to  solve  the  nonlinear  equations.  Total  time  = .194,  total 
f-evaluations  170,  total  - evaluations  = 13  7. 

' oy 

At  a moderate  tolerance  like  10  , one  correction  step  is  quite  common;  on 

harder  problems,  more  corrections  occur.  In  Table  4.5  we  illustrate  the  behavior 
of  this  code  on  a somewhat  easier  problem  (#3  of  [Daniel-Martin  (1975)]): 

ii  3 sin2rrx|,  2.  2,_  . _ . 2 sin  2 ttx  r „ , 

y y (-  y i e |4tt  (cos  2rrx  - sin  2 rrx)  - e -1J  for  0 < x ^ 1, 

/rw  /i\  , . sin  2// x 

y(0)  y (1)  1 with  solution  y(x)  e 


•rial 


TABLE  4.5  PROBLEM  #3  WITH  TOL  10  SO  THAT  THE  MAXIMUM  ABSOLUTE  ERROR  MUST  NOT 

EXCEED  2. 7 x 10~b 


, 8.3  (?  -2  . 2.8  @-9 

(2.0  <a  -2) 

j NEWT  6 

! (No  def.  corr.  step 
j possible) 


9. 6(3-4  +-  2. 9(3-9 
(1.1  & -3) 

NEWT  3 

(Not  at  correct  order) 


6.4  @ -5  + 9.4  (3  -7 
(6.4  (3-5) 

NEWT  = 2 


2.1  @ -7  + 4.5  (?  -7 
(7.1  <a  -7) 

NEWT  = 2 

(Convergence!)  STOP 


See  Table  4.4  toi  \ e .ntrics.  Total  time  .210,  total  f-evaluations 

202,  total  - eval  iati  Lt>^. 

•v 


Tables  8.  and  8.4  of  [Daniel-Martin  (1977)]  give  the  (very  similar)  be- 
havior of  NUMIDC  on  these  problems  at  these  tolerances.  We  feel  it  is  more  inter 
esting  to  consider  problems  with  y'  present  for  more  difficult  tolerances.  The 
next  example  problem  (#3YP  of  Section  7)  is  quite  badly  scaled. 


y"  lly'  + 12y  - 22e 

with  y (0)  1 and  y(15)  e15; 


the  exact  solution  is  y(x)  e so  that 


e15  3.26  x 106. 


Table  4.6  shows  the  path  taken  by  SPLIDC  on  this  problem  with  TOL  10 


TABLE  4.6  PROBLEM  #3YP  WITH  TOL  = 10  , SO  THAT  THE  MAXIMUM  ABSOLUTE  ERROR  MUST  NOT  EXCEED  3 . 26 x 10 


15 


- 12 

rhu  final  example  has  a tolerance  of  10  . As  can  be  noted  from  the  per- 

formance results,  this  is  a very  tight  tolerance  for  this  program;  we  will  give  an 
example  in  which  it  was  attained.  The  problem  we  will  consider  is  #5YP  of  Section  7; 

y"  2yy 1 

y(0)  o,  y (1)  - tan(l) 

with  exact  solution  y tan  x,  (jyjj  tan(l)  ~ 1.5. 

The  results  for  this  problem  with  TOL  10  ^ are  given  in  Table  4.7. 

This  concludes  our  brief  presentation  of  performance  data  for  SPLIDC;  more 
detailed  results  may  be  found  in  Section  7 of  [Martin- Daniel  (1977)]. 

5 . Conclusions 

The  first  conclusion  is  that  SPLIDC  does  work.  The  test  cases  mentioned 
in  Section  4 illustrate  that  iterated  deferred  corrections  with  collocation  is  a 
feasible  solution  technique  for  a substantial  class  of  problems;  however,  the 
overall  performance  of  SPLIDC  leaves  open  the  question  of  whether  collocation  methods 
with  acceleration  can  be  made  competitive  with  other  solution  techniques.  If  collo- 
cation is  to  be  made  competitive  with  the  other  techniques,  it  seems  clear  that  the 
appropriate  spline  space  must  be  carefully  chosen. 


- 


total  ^ — evaluations  194  = total  ■*—  -evaluations. 


i; 
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7 . Detailed  Numerical  Results 

In  this  section  wc  will  list  our  tost  problems  and  the  results  we  obtained 
for  them  using  SPLIDC.  We  will  list  first,  for  the  reader's  convenience,  the 
17  test  problems  not  involving  y'  from  [Daniel-Martin  (1975)).  The  sources  of 
these  problems  and  some  motivation  for  choosing  them  can  be  found  there.  The 
13  problems  used  for  Tables  4.1  and  4.2  are  those  below  numbered  1 through  12 
and  17. 


3 2 

1.  y"  y (sin  x)  (1  t sin  x)  for  0 < x < tt, 

y (0)  y (77)  0. 

Solution:  y(x)  sin  x. 

2.  y"  e^  for  0 < x < 1 , 

y(0)  y (l)  o. 


Solution:  y (x)  - In  2 t 2 In  (c  sec  (x  - ~)j 

where  c « 1.330055694906108... 


3 sin  2ttx  , 2 2 . 2 sin  2ttx 

3.  y"  y t y -i  e [4ir  (cos  2rrx  - sm  2ttx)  - e - 1J 


for  0 < x < 1, 
y(0)  y(l)  1. 


sin  2tix 


Solution:  y(x) 


20 


4.  y"  ^ ( y • x • l ) tor  0 < x v 1, 


y (0)  y(l)  0. 


Solution:  y(x)  2(2-x)  - x -1. 


’ . y”  - .000 iy/ ( .000 1 i x") “ for  - . 1 < x < . 1 , 


•y (- . l)  y(.U  .1//.0101 


Solution:  y (x)  x/\/.0001  4 x 


x/  /o 


, ,,  ! sin  4?rxfl  2 2 , , . 2 sin  4t/x  , , 

0.  y y t y i e [16rr  (cos  4irx  - sin4ux)  -e  - 1] 


for  0 < x < 1, 


y(0)  y CD  1. 


Solution:  y(x)  e 


sin  4nx 


„ 5 sin  2irx  , , 2 2 . - . 4sin  277-x  , 

/.  y y • y ■ o l 4tt  (cos  2n rx  - sin  2t?x)  - e -lj. 


for  0 < x < 1, 


y (0)  y ( li  1 


Solution:  y(x) 


o ii  ^ Sin  +ITX  * 2 2 A c { n ^7rv 

y y y ' 0 [Hat  (cos  4?rx  - sin4rcx)  - e - lj 


for  0 v x < 1. 


y(0)  y(l)  1, 


Solution:  y(x) 


sin  4 Tlx 
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3 x 2 2.  2x 

9.  y"  y i y -t  e {4f  cos  2ttx  t sin  27rx  [-4tt  - (sin  2t rx)e  ]} 


for  0 < x < 1, 


y (0)  y(i)  0. 


Solution:  y(x)  e sin  2ttx. 


2 2 

y"  400  (y  * cos  vrx)  t 2rr  cos  2?rx  for  0 < x < 1, 


y(0)  y (1)  0 


t-20x  - 20  (1-x)  2 

Solution:  y(x)  ^ cos  fx. 

1 1 e 


0 if  x < 0 or  y < 0 


11.  y"  <* 


. -b  , , 2 , 

4x  y - 6y  In  y otherwise 

.. 


for  0 < x < 1, 


y (0)  0.  y (1)  e . 


Solution:  y(x)  e 


-1/x2. 


12.  y"  1600y  - (1600  t-TT*")  sin  Trx  for  0 < x < 1, 
y(0)  y (1)  0. 


Solution:  y(x)  • sin  me. 


13.  y"  -(it  + .000 1)  y i .0001  sin  irx  for  0 < x < 1, 
y(0)  y (1)  0. 


Solution:  y(x)  sin  ?rx. 


14.  y 


-(ir“  .0000 l)y  • .00001  sin  ttx  for  0 < x < 1, 

y(0)  y (l)  0. 

Solution:  y(x)  sin  ttx. 


„ 3 40  / 1 \2 / J ( 1\8  . _ . . 

y y ' T (x_  2/  ' |^x~  2)  for  0 V X < 1, 

m8/3 

y(0)  y ( l)  . 

(1  \ 8 / 3 
X ■ "9  I 

„ 3 238  f 1 \ 1 1 / 3 ( 1 \ 1 7 _ . 

y"  y • — (x--l  '(X'2j  for  0 < X < 1, 

,(°)  yW^.p. 

Solution:  y(x) 


22  v 

17.  y"  400  (y  + cos  ttx)  + 2rr  cos  2rrx  + 400(e  -eJ  ')  for  0 v x < 1, 
y (0)  y (1)  0. 

- 20x  - 20 ( 1-x'  0 

where  y * (x)  ^ cos~Tnx. 

lie 

Solution:  y(x)  y*(x). 


We  will  now  discuss  in  more  detail  the  test  problems  we  used  that  do  involve 
y'  in  the  differential  equation. 


23 


Problem  3 1VP 

y"  for  0 < x < 2, 

l^x 

y (0)  1.  y (2)  .2 

Solution:  y 1/ (lix~) 


This  problem  was  taken  from  [pyfe  (1969)]. 


Problem  #2YP 

y"  y * -^j(y ')  esin  ” x [cos~2rrx-sin  2rrx  J- 

1 „ 3 , „ . 3 2 sin  2rrx  , •> 

■j-q  8rr  (cos  2rrx)  e - li  for  0 < x v 1 

y(0>  yd)  l. 

„ , _ . . . sin  27IX 

Solution:  y (x)  e 

This  problem  was  constructed  from  #3  by  first  substituting  y'  for  y in  the  non- 
linearity. Some  form  of  continuation  method  would  be  needed  to  handle  that  problem 
directly;  to  avoid  complications  not  relevant  to  the  present  work,  we  used  the 
factor  ^ so  that  the  present  algorithm  can  solve  the  initial  system. 

The  remaining  problems  are  from  [Scott  (1975)].  They  have  a variety  of 
interesting  characteristics. 


Problem  # 3YP 


y” 

lly  ’ 

t 12y  - 22e 

for  0 < x < 15, 

y(0) 

1, 

y(15)  el\ 

Solution: 

y(x)  eX. 

This  linear  problem  is  very  badly  scaled. 


m 


J 
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Problem  4 VP 


T i 17 

— — - ‘ 2 • tan  x f y 1 - Try 
L tan  x J j.0 


, t r , 
for  r <• 
0 


7T 

3 


No  exact  solution  known. 

The  final  three  problems  are  successively  harder  versions  of  one  basic 

problem. 

Problems  ->'5YP-7VP 

y"  2yy'  tor  0 < x v a ^ , 
y (0)  0,  y(a)  tan  (a)  . 

Solution  y (x)  tan(x) . 

As  a approaches  Tr/2,  this  problem  becomes  extremely  difficult,  since  we  are 

trying  to  follow  the  function's  behavior  as  it  approaches  a singularity.  For 

# 5YP  a 1 and  the  problem  is  relatively  simple.  For  #6YP  a 1.5  and  the 

computation  of  an  initial  solution  is  quite  difficult.  To  discover  how  the  code 

would  behave  if  it  got  started,  we  used  the  exact  solution  of  the  differential 

equation  at  the  initial  grid  points  as  an  initial  guess.  For  #7YP  a 1.55 

and  this  technique  failed  to  produce  an  adequate  initial  solution  for  any  tolerance. 

For  these  problems,  we  will  now  summarize  the  performance  of  SPLIDC 

• ^ - b 

at  each  tolerance  tested.  For  each  problem  and  each  tolerance  (TOL  10  . 10 

-9  - 12 

10  , 10  ) we  list  in  Table  7.1  through.  7.24  the  actual  errors  achieved  and 


the  cost  to  achieve  them.  (Recall  that  a requested  tolerance  of  TOL  means  that 
the  goal  is  to  make  the  error  between  the  computed  solution  Y and  the  true  solution 


2b 


r 

y*  satisfy  Y-y*  / y-  TOL.)  As  usual,  an  entry  a (3  b denotes  a«10“. 

If  SPLIDC  detects  that  the  requested  error  cannot  be  attained,  we  will  list 
the  condition  detected  as  well  as  the  error  we  were  able  to  attain.  (The  listed 
function  counts  include  the  unnecessary  evaluation  at  the  endpoints,  as  dis~ 
cussed  in  Section  4.) 


TABLE  7.3  Performance  Data  for  Problem  #3 


TABLE  7.4  Performance  Data  for  Problem  #4 


TABLE  7.5  Performance  Data  for  Problem  #5.  y*  - .995 


TABLE  7.6  Performance  Data  for  Problem  3A6  y*  - e = 2.718 


TABLE  7.9  Performance  Data  for  Problem  #9  y*  « 2.144 


TABLE  7.10  Performance  Data  for  Problem  ?*10 


TABLE  7.11  Performance  Data  for  Problem  // 1 1 . i;y*;  e ~ .168 


TABLE  7.13  Performance  Data  for  Problem  #13 
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TABLE  7.15  Performance  Data  for  Problem  #15 
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TABLE  7.18  Performance  Data  for  Problem  #1YP 


TABLE  7.19  Performance  Data  for  Problem  #2YP 
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TABLE  7.20  Performance  Data  for  Problem  #3 YP 


TABLE  7.22  Performance  Data  for  Problem  5YP 


TABLE  7.23  Performance  Data  for  Problem  #6YP.  jiy*|;  tan  15 


TABLE  7 . 24  Performance  Data  for  Problem  //7YP.  y*  = tan  1.55  » 48. 
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2B  Ji-L  77 

SPLIDC  - 1 < 

ODDD  ’ CC 

DO  C C 
D D C 

D 0 C 

D D C 

D D C r 

ODDD  CCC 


PUMPOUT  INF  SPIIDC<F,DFY,DFYP,N.A.H,AlFO.BVO,GO,AlFN,BVN.GN.W,X, 

] TOL .MO YP . Y I NT ) 

C THIS  PROGRAM  IMPlFMf-NTS  ITFPATFD  DEFERRED  CORRECTION 

C ON  A SPlinE  collocation  method  FOP 
C 

C (O'**?)  Y^F  ( X . Y.YP)  » 

C 

c with  genfpal  sepaparlf  linfar  boundary  conditions. 

c 

C ALFO*Y (A) «BVO*YP ( A) =G0  . ALFN*Y ( B) .fiVN* YP ( B) =GN 

C 

C AS  DESCRIBED  IN  A FORTHCOMING  REPORT. 

C 

C While  ThE  METHOD  SUPPORTS  GENERAL  LINEAR  BOUNDARY  CONDITIONS 

ThE  PRESFNT  IMPLEMENTATION  DOES  NOT  SUPPORT  BOUNDARY  CONDITIONS 
C INVOl  VlNG  YP. 

C 000000*00  EXTRANEOUS  FEATURES  ooooooooooo 

THIS  code  CONTAINS  A NUMBER  of  features  oriented  toward 
RESEARCH  WORK  rather  than  production  work,  such  as  trace-type 
PRINTING  AND  ESTIMATION  OF  COMPUTATIONAL  EFFORT.  THE  VARIABLE 
C IFLAG  IS  DECLARED  AS  COMMON  AND  TRACE-TYPE  PRINTING  OCCURS 

IF  AND  ONLY  IF  IFLAG  = 0.  T HE  COMMON  VARIABLES  K0UNT1  AND  K OUNT 2 
COUNT  T Hf  NUMBER  OF  FUNCTION  EVAilJATIONS  PERFORMED. 

C SPECIAt  TRACE  PRINTING  IS  DONE 

C IN  THE  ROUTINE  PRSUB  . ALL  OF  THESE  ITEMS  AND  REFERENCES 

C TO  Them  CAN  BE  REMOVED  IN  A PRODUCTION  CODE. 

C 

C oooooo  DESCRIPTION  OF  SUBROUTINE  PARAMETERS  oooooooo 

c 

C F : AN  EXTERNAL  USER-PROVIDED  FUNCTION  (OF  TYPF  REAL) 

C WITH  CALLING  SEQUENCE  F(X,Y,yP)J  THIS  SHOULD  COMPUTE 

C The  FUNCTION  F DEFINING  THE  DIFFERENTIAL.  EQUATION 

C FOR  SIMPLE  VARIABLES  X.Y.YP. 

C OF  Y : AN  EXTFRNAL  USER-PROVIDED  FUNCTION  (ALSO  OF  TYPE 

C REAL)  WITH  CALLING  SEQUENCE  DFY(X.Y.YP)  WHICH  COMPUTES 

C THF  CORRESPONDING  VALUE  OF  OF/DY. 

C OF  YP  : ANOTHER  EXTFRNAL  USER-PROVIDED  FUNCTION  (ALSO  OF  TYPE 

C REAL)  WITH  CALLING  SEQUENCE  DFYP(X.Y.YP)  WHICH  COMPUTES 

C THE  CORRESPONDING  VALUES  OF  DF/DYP. 

c if  noyp  is  true,  dfyp  IS  never  CALLED. 

C N : AN  OUTPUT  PARAMETER.  THE  FINAL  APPROXIMATE  SOLUTION  IS 

C PROVIDED  ON  A UNIFORMLY  SPACFD  MESH  OF  WIDTH  H=(XF-XO)/N. 

C AS  PRESENTLY  WRITTFN  THIS  CODE  IS  LIMITFD  TO  N LESS  OR 

C EQUAL  P56J  THIS  CAN  RF  CHANGEO  By  APPROPRIATE  CHANGES 

C IN  DIMENSION  STATEMENTS  ANO  IN  THE  VALUES  OF  NMAX(=?56) 

C AND  NMAXH ( = 12H) . 
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It-  VINT  IS  TRUE.  N IS  ALSO  AN  INPUT  PARAMETER. 

A : LEFT  END  POINT  OF  INTEGRATION  I NTER V AL ♦ PROV I OE D BY  USER. 

B : RIGHT  END  POINT  OF  INTEGRATION  INTERVAL.  PROVIDED  BY 

user;  a = b is  an  invalid  input. 

ALF0.BV0.G0  ' COEFFICIENTS  OF  BOUNDARY  CONDITION  EQUATION  AT  A. 

ALFN.BVN.GN  - COEFFICIENTS  OF  BOUNDARY  CONDITION  EQUATION  AT  9. 

RVO  AND  BVN  BOTH  ZERO  REMOVE  YP  EROM  THE  BOUNDARY  CONDITIONS. 

IT  IS  THIS  CASE  WHICH  IS  FULLY  IMPLEMENTED.  THE  PRESENT  THEORY 
REQUIRES  THAT  THE  LINEAR  INTERPOLATION  PROBLEM  HAVE  A UNIQUE 
SOLUTIONS  IE  THIS  IS  NOT  TRUE.  EXECUTION  IS  TERMINATED. 

W : AN  OUTPUT  ARRAY,  CONTAINING  THE  COMPUTED  APPROXIMATE 

SOLUTION.  UPON  EXIT,  THE  FIRST  FOUR  COLUMNS  OF  W CORRESPOND 
TO  THE  ELEMENTS  OF  X.  THEIR  ELEMENTS  ARE  APPROXIMATIONS 
TO  THE  FUNCTION  VALUF  AND  THE  FIRST  THREE  DERIVATIVES  OF  THE 
SP|  INE  SOUGHT  BY  THIS  NUMERICAL  METHOO. 

IE  YINT  IS  TRUE.  W IS  ALSO  AN  INPUT  PARAMETER. 

X : AN  OUTPUT  ARRAY  CONTAINING  THE  N ♦ I POINTS,  INCLUDING 

A AND  B,  AT  WHICH  THE  SOLUTION  IS  APPROXIMATED.  UPON 
EXIT,  SPLIOC  WILL  HAVE  XU)  = A , X(NU)  = B.  AND 
X(I)  = A ♦ (I-I)*H  , WHERE  H = (6  - A)  / N 

TOL  : REQUESTED  ACCURACY,  PROVIDED  BY  USER.  SPLIDC  ATTEMPTS 

TO  COMPUTE  AN  APPROXIMATE  SOlUTION  Y WITH  A MAXIMUM  ERROR 

less  than  tol  times  the  maximum  computed  value  of  y,  a 

SORT  OF  SCALED  MAX-NORM  ERROR.  TOL  IS  ALSO  AN  OUTPUT 
PARAMETER,  CONTAINING  ON  EXIT  FROM  SPLIDC  THE  ESTIMATED 
UPPFP  BOUND  ON  THE  ERROR.'  THIS  ESTIMATE  IS  THE  SUM  OF  THE 
ASYMPTOTIC  ESTIMATE  OF  THE  DISCRETIZATION  ERROR  AND  THE 
Estimate  OE  THE  ERROR  IN  SOLVING  THE  NON-LINEAR  EQUATIONS. 

AN  INPUT  VALUE  OF  LESS  THAN  I.E-12  WOULD  LEAD  TO  UNRELIABLE 
PERFORMANCE  ON  THIS  MACHINE.  SO  SUCH  VALUES  ARE  REPLACED 
BY  1.E-I2  WITH  AN  APPROPRIATE  WARNING  MESSAGE  TO  THE  USER. 

THF  NUMBER  l.E-12  IS  APPROPRIATE  FOR  THE  CURRENT 
CONFIGURATION.  IN  GENERAL  WE  WOULD  RECOMMEND  A NUMBER 
APPROXIMATELY  EQUAL  TO  THE  PRODUCT  OE  THE  MAXIMUM  NUMBER  OE 
POINTS  AND  THE  PRECISION  OE  THE  MACHINE. 

NO YP  : LOGICAL  VARIABLE.  IE  NOYP  IS  TRUE,  YP  IS  NOT  INVOLVED 
IN  EITHER  THE  DIFFERENTIAL  EQUATION  OR  THE  BOUNDARY 
CONDITIONS,  SO  THE  ROUTINE  SPLNYP  IS  USED.  OTHERWISE. 

THE  ROUTINE  SPL.WYP  IS  USED. 

YINT  : IF  YINT  IS  TRUE,  THF  USER  SUPPLIES  AN  INITIAL  GUESS 

TO  The  SOLUTION  Y.  (OTHERWISE,  THE  CODE  BEGINS  WITH 
A LINEAR  INTERPOLATION  OE  THE  BOUNDARY  DATA.)  IE  YINT 
IS  TRUE  N IS  MEANINGFUL  ON  INPUT.  THE  FIRST 
N * 1 ELEMENTS  CE  THE  FIRST  COLUMN  OF  W SHOULD  CONTAIN 
THF  DESIRED  INITIAL  APPROXIMATION  TO  THE  SOLUTION. 

SINCE  THE  OIFFEcENCF  IN  THF  CODE  WITH  AND  WITHOUT 

rp  is  substantial.,  we  simply  branch  to  the  appropriate 
routine  depending  on  noyp. 

EXTERNAi  f.dey.dfyp 
DIMENSION  w <?59, S) , x <?SB) 

LOGICAL  noyp.yint 
IE  11*  0 YP  ) GOTO  ?0 

CALL  SPLWYP (E  ,DEY,DEYP.N,A,B, AlF0.BV0.G0. ALFN.RVN.GN.W.X, 

1 T QL  ♦ Y I N T ) 

RETURN 
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20  C^LL  c-PLNYP(F»DFy«DFYP.N»A»R,ALFO»BVO,GO»AlFN»BVN»GN.W.X» 

1 TOL.YINT) 

RETURN 

END 


SUBROUT  INF  SPl  NYP <F .DFY.DFYP.N.A.B.ALEO.BVO.GO.ALFN.BVN.GN.W.X. 

1 T OL , Y ’ NT ) 

wF  ALLOCATE  STORAGE  AT  THIS  POINT.  w IS  USED  TO  STORE  THE 
MATRIX  OF  THE  LINEARIZED  SYSTEM  DURING  EXECUTION. 

FU.E0.AND  El  STORE  THE  COMPUTED  VALUES  OE  F DF/DY  AND  DF/DYP,  RESP. 
THE  ARRAY  S IS  PRIMARILY  USED  TO  STORE  COMPUTED  APPROXIMATIONS 
TO  THE  DISCRETIZATION  ERROR;  THE  ARRAY  RHS  IS  USEO  IN  THE 
SOLUTION  OF  THE  NON-LINEAR  SYSTEMS  RHS  CONTAINS  THE  COMPUTED 
RESIDUAL  WHICH  IS  overwritten  WITH  THE  CHANGE  IN  Y. 

WK  IS  A WORK  AREA  FOR  The  LINEAR  EQUATION  SOLVER. 

RHSTOT  CONTAINS  THE  SPLINE  COEFFICIENTS  FOR  THE  COMPUTED 
SOLUTION.  Y.YP  AND  AY2P  CONTAIN  COMPUTED  APPROXIMATIONS 
TO  Y.YP  'ND  THE  DIFFERENCED  SECOND  DERIVATIVE  OF  THE 
COMP.  SOLUTION  RESPECT  IVELY  . 

DIMENSION  W (259,5)  . X (259)  , RhS  (259)  , E 0 < 259)  ,E1(2S9),W<(777) 

DIMENSION  RHSTOT (259)  , Y (259)  , YP (259)  . ft Y 2 P (259)  ,EU (259)  .S (2S9) 

WE  WILL  USE  THE  ARRAYS  ERRSAV.  TESTl,  AND  TES T2  IN  THE  TEST 
FOR  CORPS C T ORDER.  THE  SOLUTION  AT  LEVEL  K=0.1.2.3»4,  OR 
5 SHOULD  REhavE  at  OROER  4.8,12.  16.2').  OR  24  RESP.  WHEN 
THE  MESH-SPACING  IS  HALVED  THE  ERROR  ESTIMATE  SHOULD  DECREASE 
BY  FACTORS  OF  16. 256. 4096. FT  CETERA.  THE  ELEMENTS  OF 
TESTl  AND  TEST2  CONSTITUTE  (RATHER  COARSE)  LOWER  AND 
UPPER  BOUNDS.  RESPECTIVELY.  ON  WHAT  IS  AN  ACCEPTABLE  LEVEL 
OF  DECREASE,  WHILE  ERRSAv  SAVES  THE  ERRORS  FROM  THE  PREVIOUS 
N TO  ALLOW  COMPUTATION  OF  THESE  FACTORS  AND  ALSO  TO  HELP 
DE  TERM  I ns  THE  LEVEL  AT  WHICH  INTERPOLATED  VALUES  ARE 
ACCEPTED  WHEN  N Is  DOUBtED. 

DIMENSION  ERRSAV <6),TEST1(6),TEST2<6> 

DATA  TESTl/8.,64.,512.,4096. .3.E4.2.E5/ 

DATA  TEST?/ 32. . 1 024. ,3.F4,1.E6,3.F7,1.E9/ 

WE  WILL  USE  LOCAL  LOGICAL  VARIABLES: 

DIVNPW  (TRUE  WHEN  NEWTON^S  °E T HOD  HAS  DIVERGED) 

NEWN  (TRUE  WHEN  WE  ARf  COMPUTING  AT  A NEW  VALUE  OF  N) 

RERROR  (TRUE  WHEN  WE  HAVE  DETECTED  DOMINANT  ROUNDING  ERROR) 

CON V OB  (TRUE  WHEN  WE  HAVE  HAD  CONVERGENCE  OF  NEwTONxS  METHOD 
AT  LEAST  ONE  TIME.) 

NOEVAL  (TRUE  WHEN  E AND  DFY  HAVE  ALREADY  BEEN  EVALUATED.) 
LOGICAL  yint .rerror.di VNEW.CONVOR.NEWN.NOEVAl 
common  /FLAG/  ielag 
COMMON  /KOUNT/  KOUNT 1 .KOUNT2 
PATA  N A X /?56/ 

DATA  Nv  A x h /128/ 

CHECn  The  VALIDITY  OF  The  input  PARAMETERS 
I F ( A .EO.  B ) ST OP 2 6 

9.99F-13  EQUALS  l.E-12  - FPS  FOR  OUR  PURPOSES 
IF ( TOL  ,LF.  9.99F-1 3) STOP30 
IF  ( .NOT.  Y I N T ) N = 8 

INITIALIZE  SOME  VARIABLES 
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The  INITIAL  value  or  N is  stored  as  no;  when  n=no  NO  BACKGROUND 
INFORMATION  fxists. 

MAXr  is  THE  HIGHEST  PREVIOUS  LEVEL  OF  DEFERRED  CORRECT  ION  THAT 
HAS  BEEN  REACHED.  WHILE  K IS  THE  PRESENT  LEVEL. 

EPNEWT  IS  THE  TOLERANCE  FOR  CONVERGENCE  IN  NEWTON-*S  BE  T HOD . 

2 MAXK=0 
N0=M 

E PSN  W T = T 01  /2. 

RERROR=. FALSE . 

CON VOB= . F ALSE  . 

-=0 

BEGINNING  OF  COMPUTATION  WITH  A NEW  VAlUE  OF  N. 

3 ERR-I  =1.E20 
NP 1 =N  ♦ 1 
NP  2 = N *2 
N P 3 = N «■  3 

-MAX  IS  MAXIMUM  ALLOWABLE  NUMBER  OF  EXTRAPOLATIONS  FOR  A 
GIVEN  N.  THERE  APE  THREE  LIMITS.  THERE  MUST  BE  ENOUGH  POINTS 
TO  COMPUTE  THE  EXPANSIONS;  6 IS  DETERMINED  BY  THE  SIZE  OF  THE 

relevant  arrays;  maxk*2  is  a restriction  on  the  number  of 

EXTRAPOLATIONS  WITHOUT  HISTORY. 

KMAX=  (NP  1 -A ) /A 
KMAX=MTN0(KMAX,6.MAXK*2) 

NFWN=.TRUE . 

EO ( 1 ) r ALFO 
El ( 1 ) = RVO 
EO  <NP3) =ALFN 
El  (NP3>  =RVN 
( B — ) /N 

WF  START  PREPARING  TO  USE  NEWTON-S  METHOD  TO  SOLVE  THE 
NON-LINEAR  EQUATIONS  APPROPRIATE  TO  THE  SPLINE  COLLOCATION  METHOD 
USED.  THE  VARIOUS  CORRECTION  TERMS  MUST  BE  INITIALIZED 
APPROPRIATELY.  TO  ZERO  AT  LEVEL  ZERO. 

ITLIM  AND  ITLIM?  LIMIT  THE  NUMBER  OF  NEWTON  ITERATIONS  AS 
DESCRIBED  LATER. 

I Tl.  IM=  1 n 
I T L I M2  = 1 3 
nO  A T -2  .-•P2 
6 F 1 ( I ) = 0 . 

x ( 1 ) = A 
< (NP1 ) =B 
nO  10  1=2. N 

10  X ( I ) = x ( 1 ) ♦ ( I - 1 ) *h 
00  12  1=1. NP3 
12  S U)  = 0 . 

IF(N  ,rQ.  NOIGOTO  1 R 
C FILL  IN  Y values  after  DOUBLING  N 

CALL  GFTVAL(N.H,RHFT0T.Y,YP.AY?P) 

c after  an  interpolation,  we  use  the  previously  computed  values 

c OF  F and  DF/DY.  the  next  few  lines  HANDLE  that. 

NHAI  f =N/2 
N I =NH  A(  F ♦ 1 
N2  = NHA|  F ♦ ? 

FULNPl ) =FU ( N 1 ) 

FO (NP2) =E  0 ( N2 ) 

DO  1 A 1 = 1 .NHALF 
NO  I =N 1 - I 
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NI =NP 1 -2<>  I 
^ U ( N I ) =FU (NOI  ) 

N I =N I ♦ 1 

EO (Ml ) = E0 (NOI  ♦ 1 ) 
p 1 1 ( N I > = F ( X ( N I ) t Y ( N I ) ♦ Y P ( N I ) ) 

14  E 0 ( N I ♦ 1 ) =-DF  Y ( X ( N I ) tY (NI)  * YP ( N I ) ) 

K OUNT 1 =KOUNT 1 *NHALF 
KOUNT?=KOUNT2*NHALF 
I F ( K ,PQ.  0 ) GOT  0 34 

CALL  maTMAN(W,RHS,FO,E i ,WK,N,H, 1 ) 

IF  ( I FI  AG  .GT.  OIGOTO  1002 
CALL  PPSUR <N»28fERRN0R»X,RHST0T ) 

1002  CONTINUE 

NEXT,  EVALUATE  THE  NEW  CORRECTION  4 FOR  THE  NON-LINEAR 

equations  that  must  re  solved  at  the  new  level  k after  an 

INTERPOLATION  FOR  THE  VALUES  ACCEPTED  AT  LEVEL  K-l. 

WHEN  NFWN  IS  TRUE  , DEY  IS  EVA'UATED  DURING  THE  NEWTON 
ITERATIONS  . 

NEWN= .FALSE . 

I TL IM- 3 
I T |_  iM-’-S 

CALL  SPLDCG(x ,NtFU»S, IERPOR) 

HO  1 R I=i,K 

18  EWRSAV  ( I ) =FRRSAV  ( I ) “16.°"°  (-1  ) 

OQT0  <8 

10  I F ( .NOT.  Y I NT ) C ALL  L I NST ( GO • AL  FO  « GN  * ALEN , A » B * N ♦ H » RHSTOT ) 

IF ( Y I NT ) C ALL  GETST ( W . RHSTOT . WK , N , H ) 

I F ( I FI  AG  .EQ.  0 ) C ALL  PRSUB(N,?4,TW0PI.X, RHSTOT) 

C M A I N PORTION  OF  NEWTON-*S  METHOD 

C WF  ENTER  here  WITH  an  initial  guess  for  Y,  EITHER  the  LINEAR 

C INTERPOLATION  OF  THE  BOUNDARY  DATA  OR  AN  IMPROVED  ESTIMATE 

C FROM  a LOWER  LEVEL  AT  THIS  N OR  A HIGHER  LEVEL  AT  A SMALLER  N 

C OR  A USER  SUPPLIED  INITIAL  GUESS. 

C AT  LEVEL  k = 0 WE  ALLOW  11  ITERATIONS  BEFORE  TERMINATING  THE 

C PROCEDURE.  AT  HIGHER  LEVELS,  ONLY  4 ITERATIONS  ARE  ALLOWED. 

C AT  LFVFL  k-0  NEWN  IS  TRUE  AND  WE  EVALUATE  DFY  AT  EACH 

C ITERATE.  IN  MOST  CASES  F AND  DFY  ARE  EVALUATED  BEFORE  ENTERING  THE 

C LOOP  AND  ARE  NOT  REEVALUATED  ON  THE  FIRST  ITERATION?  IF  NO  SYSTEM 

C IS  SOLVED  AT  K=0  THIS  IS  THE  ONLY  TIME  DFY  is  EVALUATED. 

C NORM/-  EXIT  OF  THF  NEWTON  ITERATION  OCCURS  WHEN  THE  NORM 

C (PNORM)  OF  THE  CHANGE  IN  THE  SOLUTION  Y IS  LESS  THAN  OR 

C E QU A(  TO  HALF  THE  USER^S  TOLERANCE  TIMES  THE  NORM  OF  Y. 

C abnormal  EXITS: 

C (1)  THE  NORM  OF  THE  RESIDUAL  INCREASES  FROM  STEP  N TO 

C N ♦ 1,  WHERE  N > 0. 

C nR 

C (2)  THE  NUMBER  OF  ITERATIONS  EXCEEDS  ITLlM  AND  CONVERGENCE 

C IS  NOT  EXPECTED  RY  ITLIm>. 

C IF  THF  CHANGE  IN  Y (RNORM)  INCREASES  BUT  THE  RESIDUAL  DECREASES, 

C THIS  IS  TAKEN  AS  A SIGN  THAT  ROUNDING  ERROR  IS  PREVENTING 

C CONVERGENCE  TO  THE  TOLERANCE  AND  THE  PROGRAM  GIVES  AN  ERROR 

C MESSAGE  AND  TERMINATES. 

38  1 N n E X = 0 

DIVNEW-.f  Al  SE  . 

NnFVAL  = . TRUE . 

I F < ( N . EQ . NO)  .AND.  (k  .EQ.  0)  )NOEVAL  = . false. 

40  CALL  GE TVAL <N»H. RHSTOT, Y,YP,AY2P) 
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Phs ( 1 ) =0 . 0 
R^S  (*  P 3 ) =0 . 0 

RESN=A~  AX  1 ( ABS  (WHS  ( 1 ) ) . ABS  (WHS  (MP3  ) ) ) 

IF  (NOEVAI.  ) GOTO  6B 
HO  BO  I = 1 ,NP 1 
YP  TF  MP  = 0 . 0 

FU( I > =F (X ( I > »Y  < I ) ,YPTE  MP) 

I F < .MOT.  NLWN)GOT0  60 

EO  < I ♦ 1 ) =-DFY  (X  < I ) , Y( I ) , YPTEMP) 

T0Mr IMUF 

KOUMTl=KOUMT 1 ♦ MP 1 
IF ( .MOT . MFWM) GOTO  66 
KOUNT  2=K0UNT2*NP1 
NOEVAL=. FALSE. 

60  70  1=1. MP1 

RHS( I ♦ 1 ) =S ( I ♦ l ) »FU( I ) -AY2P ( I» 

IFfPFSM  .LT.  ABS (RHS < I ♦ 1 > ) ) RESN=ABS (RHS ( I ♦ 1 ) ) 

CONTINUE 
T JO«  = ? 

IF (NEWN) I JOB=0 

CALL  MATMAN (W.RHS.E0.E1 .WK.N.H. I JOB) 

RNORM=0 . 

RTNORM=0 . 

no  100  1=1, MP3 

RHSTOT ( I ) =RHSTOT ( I ) ♦ RHS ( I ) 

I F ( RNORM  .IT.  ABS (RHS ( I ) ) ) RMORM=ABS (RHS ( I ) ) 

IF(«TNORM  .LT.  ABS (RHSTOT ( I ))) RTMORM=ABS (RhSTOT ( I > > 

AYNORM  IS  AN  APPROX  I M A T I ON  TO  T HF  NORM  OF  Y BASEO  ON  THF 
THEORETICAL  BOUND  ON  THE  RATIO  OF  THE  NORMS  OF  THE  SPLINE  COEF 
AND  THE  FUNCTION  Y. 

A YNORM=RTNORM/3. 

INDFX= INDEX ♦ 1 

I F ( I Nn£  x . EO . l)RNORMl=RNOR*< 

FPS=EPSNWT*AYNOR^ 

IF(RN0RM  .LE.  EPS)G0T0  3Q0 

IF ( INDEX  .LE • 2)  GOTO  120 

IF ( RE SN  .GT.  RSNSAV1GOTO  300 

IF((RN0Rm  .GT.  RN6AV)  .AND.  CONVOBIGOTO  350 

RNS?=RNSAV 

PNS AV=RN0Rm 

RSNSAV=RESN 

I F ( I NDE  X .LF.  ITLIMJGOTO  AO 
IFIINDEX  .GT.  I TL  I M<? ) GO T 0 300 

I F ( ( RNORM  « (RNQRM/RNS2)  * * ( I TL  I M<f  - I N EX))  .LT.  EPS)G0T0  AO 
NFWTON-S  METHOD  IS  DIVERGING  AND  WE  MUST  DOUBLE  N IF 
POSSIBLE.  IF  THE  FINAL  NEWTON  ITERATE  SOLVES  THE  NON- 
l INFAR  SYSTEM  BETTER  THAN  THE  FIRST  ITERATE.  THEN  WE  GET 
AN  ERROR  ESTIMATE  FO»  USE  IN  TESTING  ORDER  AND  WE  USE 
V A | UE S INTERPOLATED  FROM  ThF  LAST  ITERATE  IN  ORDER  TO 
START  NFWTON-S  METHOD  AT  LEVEL  /FRO  AND  WITH  N DOUBLED  J 
IF  THE  FINAi  ITERATE  WAS  WORSE.  HOWE  VF  R » WE  START  WITH 
ThF  LINFAR  INTERPOLATION  OF  THE  BOUNDARY  DATA. 

D I VNE  W= . TRUE . 

I r ( ( RNORM 1 .GT.  RNORm)  .OR.  CONVOH)GOTO  A00 
I F ( N .GT.  NMAXH)GOTO  BOO 
IF (YlNT)GOTO  ft? 0 
N — 
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3B0 


400 


450 


4*0 


100  1 
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goto  ? 

WF  NOW  THINK  THAT  ROUNDING  ERROR  HAS  CONTAMINATED  OUR 
PFSULTS.  SO  Wf  PREPARE  TO  QUIT  COMPLETELY  AFTER  ESTIMATING 
THE  ERROR  WF  DID  ACHIEVE. 

RfRROR- . TRUF . 

GOTO  400 

WE  NOTE  FOR  FUTURE  REFERENCE  THAT  CONVERGENCE  HAS  BEEN 
OBTAINED  AT  LEA  ^T  ONCE. 

CO  V08  = .TRUE » 

WF  WANT  TO  ESTIMATE  ThE  D I SCRE T I Z A T ! ON  ERROR » ESSENTIALLY 
BY  ONE  STEP  OF  NEWTONr*S  METHOD.  WE  SET  UP  THE  RIGHT  HAND 
SIDE  RHS  NEEDED  FOR  THE  LINEAR  EQUATIONS  OF  THIS  ONE  STEP  5 
WHILE  COMPUTING  RHS  WE  WILL  ALSO  COMPUTE  S.  THE  CORRECTION 
(ESTIMATING  LOCAL  TRUNCATION  ERROR)  THAT  WE  WILL 
NPED  FOR  THE  NON-LINEAR  EQUATIONS  AT  THE  NEXT  LEVEL  UP  IF 
WE  GET  THAT  FAR. 

K =K  ♦ 1 
T T I I M = 3 
I T'  I M2=-5 

Call  SPLDCG(K,N,FU.RHS» I ERROR) 
p • s ( 1 ) r o . 0 
R- ' (MP3) =0.0 
DO  450  1=2. NP2 
S T F =RHS ( I ) 

RHS ( I ) =S ( 1 ) -STF 
S ( I ) =STF 

WITH  ThE  NEW  RIGHT  HAND  SIDE.  WE  GENERATE  AND  SOLVE  THE 
NEW  SYSTEM. 

CALL  MATMANIW.RHS.EO.El , WK ,N»H.?) 

ERROR  CONTROL  ANO  DECISION  CENTER 

Rhf  NOW  CONTAINS  OUR  ESTIMATE  OF  DISCRETIZATION  ERROR. 

WE  COMPUTE  ITS  NORM  IN  ERRNOR.  ALLOWING  US  TO  TEST  FOR 
CONVERGENCE  OE  THE  OVERALL  ALGORITHM. 

NEWN=. FALSE. 
f OR-  j0R  = 0 . 

V NORMrO . 

r,0  460  1 = 1.  IP  1 

T F = A8S ( ( RHS ( I ) *4.<RHS (1*1)  .RHS  < I *2) )/6. ) 

I F ( T E .GT.  ERRNOR) ERRNOR  = TF 
I F ( VNORM  ,lT.  ARS ( Y ( I ) ) ) YNORM=ABS ( Y ( I ) ) 

IF  ( I F L A G .GT.  0)GOTO  1001 
K 1 :K-  1 

PRINT  470. INDEX. RE SN.RNORm 
CALL  RRSUB(N,k1 . ERRNOR. X.RhS TOT) 

COMT INUF 

next  WF  TEST  OUR  ESTIMATE  OE  OVERALL  ERROR  AGAINST 

The  USF.R^S  TOLERANCE.  ERRNOR  ESTIMATES  The 

DISCRETIZATION  ERROR  FOR  ThF  EXACT  SOLUTION  OF  THE 

NON-LINEAR  SYSTEM.  WHILE  RNORM  ESTIMATES  THE  ERROR  IN 

SOLVING  THAT  NON-LINEAR  SYSTEM.  ERREST  ESTIMATES  THE 

TOTAl  FRROR  undfr  the  pessimistic  ASSUMPTION  THAT 

ThF  OTHER  ERRORS  SUM. 

eppfst=eprnor*rnorm 

TOLLOC  = TOLuYNORm 

TEFT  For  SUCCESS 

IE  'ERREST  .LT.  T 0LL OF ) GOTO  OOO 


A 
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C TERMINATE  ALGORITHM  IF  RERROR  TRUE 

IF  f Rr  RRQR ) GO T 0 9B0 

I F K .LF.  MAXK  (GO  THAT  THERE  IS  HISTORY  TO  COMPARE 
AGAINST)  WE  NEXT  CHECK  WHETHER  THE  ERROR  IS  BFHAVING 
AT  THE  CORRECT  ORDER. 

If (K  .GT.  MAXK ) GOTO  4H0 
RAT=ERRSAV(K) /ERRNOR 
FPRSAV (K ) =ERRNOR 
IF ( I F L AG  .EO.  0 ) PR  I N T 475. RAT 
This  IS  THE  TEST  FOR  CORRECT  ORDER. 

IE(RAT  .LT.  TESTl(K)  .OR.  RAT  .GT.  TEST2 <K ) ) GOTO  700 
480  ERRs AV (K ) =errnor 

The  NEXT  STATEMENT  TESTS  WHETHER  THE  PREVIOUS  DEFERRED 
CORRECTION  IMPROVED  THE  ACCURACY  ACCEPTABLY  AND 
WHETHER  AN  ERROR  ESTIMATE  WOULD  BE  ALLOWED  AFTER 
The  next  CORRECTION. 

I E ( ERRNOR  .GE.  .l*ERROLD  .OR.  K ♦ 1 .GT.  KMAX ) GOT  0 700 
CHECK  for  NEWTON  DIVERGENCE 
IF (DIVHEW)GOTO  700 
ERROLD=ERRNOR 

GO  AHEAD  AND  SOLVE  THE  SYSTEM  AT  THE  HIGHER  ORDER 

GO to  38 

WE  NEED  TO  INCREASE  N AND  THE  NEXT  SECTION  OF  CODE  HANDLES 
THAT.  WE  FIRST  COMPUTE  THE  LEVEL  INTERPOLATED  VALUES  WILL  BE 
ACCEPTED  AT.  THEN  WE  DO  THE  INTERPOLATION. 

700  I F ( N .GT.  NMAXH)GOTO  800 
M A tK 
NOL'=N 
M r ?<>N 
K TN  r -K 

ERROl.D- AM  1 N 1 (ERROLD. ERRNOR) 
pECIDF  LEVEL  TO  INTERPOLATE  TO 
00  710  1 = 1 .K 

IF(ERROLD  .LT.  ERRSAV ( I ) (-1 ) ) GOTO  710 
GOTO  7?0 
710  CONTINUE 

720  K = I - I 

IF (01  NEW) K=0 

INSPL  PERFORMS  THE  INTERPOLATION  AND  FILLS  IN  Y VALUES. 

CALL  INSPL  (K INT  »N*Y .W.S.RHSTOT  « WK  * K » H ) 

.GOTO  T 

this  is  the  n too  rig  exit 

800  PR  I ' 1 T 8 1 5 

GOTO  900 

820  IS  THE  ERROR  EXIT  TAKEN  WHEN  AN  INITIAL  GUESS  FAILS. 

820  PRINT  825 

THE  FO (LOWING  CODE  SEGMFNT  PREPARFS  TO  RETURN  TO  THE  USER. 

IN  THE  ABSENCE  OF  EPROR  MESSAGES.  THE  RETURN  IS  SUCCESSFUL. 

900  T OL  =FRRF  c T 

IF ( I F L A G .FO.  0 ) PR  I N T 920.YNORM 
CALL  SPLGFN (W.N.H.RHSTOT ) 

return 

C 980  IS  THE  ROUNDING  ERROR  EXIT. 

980  PRINT  990 
GOTO  900 

A 70  FORMAT (I10.2F20.10) 

A 75  FORMAT (80X  .F 12.4) 
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815  rcRMAr<» 
«25  FORMAT  ( » 
920  FORwAT(« 
990  f OP'"'  A T ( •» 
END 


N TOO  BIG*) 

illegal  to  reinitialize 

YNORM  = » . f 20 . 1 0 ) 
ROUNDING  ERPOR  EXIT®) 


WHEN 


INITIAL  guess  GIVEN®) 


SUBPOUT  INE  SPLWYP (F .DrY.DFYP.N. A.R. ALF0.BV0.G0. ALFN.BVN.GN.W.X. 

1 T OL  » Y t mT ) 

C WE  ALLOCATE  STOPAGF  AT  THIS  POINT.  W IS  USEO  TO  STORE  THE 

C MATRIX  OF  Thy  LINEARIZED  SYSTEM  DURING  EXECUTION. 

C FU.EO.AMD  Ei  STORE  THE  COMPUTED  VALUES  OF  F DF/DY  AND  DF/DYP,  RFSP. 

C THE  ARRAY  S IS  PRIMARILY  USED  TO  STORE  COMPUTED  APPROXIMATIONS 

C TO  THE  DISCRETIZATION  ERROR?  THE  ARRAY  RHS  IS  USED  IN  THE 

C SOLUTION  OF  THE  NON-LINEAR  SYSTEM?  RHS  CONTAINS  THE  COMPUTED 

C RESIDUAL  WHICH  IS  OVERWRITTEN  WITH  THE  CHANGE  IN  Y. 

C WK  IS  A WORK  AREA  FOR  THE  LINEAR  EQUATION  SOLVER. 

C RHSTOT  CONTAINS  THE  SPLINE  COEFFICIENTS  FOR  THE  COMPUTED 

C SOLUTION.  Y.YP  AND  A Y2P  CONTAIN  COMPUTED  APPROXIMATIONS 

C TO  Y.YP  amD  the  differenced  second  derivative  oe  the 

C COMP.  SOLUTION  RESPECTIVELY. 

DIMENSION  W (259,5) . X (259) ,RHS (259) ,E0 (259) »E 1 ( 259 > . W« ( 777 ) 

DIMENSION  RHSTOT (259) , Y (259) , YP ( 259) , AY2P (259) .EU (259) .S (259) 

C YPC  CONTAINS  THE  CURRENT  VALUES  OF  THE  YP  CORRECTION.  YPCOLD 

C CONTAINS  THE  PREVIOUS  VALUES. 

DIMENSION  YPC (259) ,yPC0LD(259) 

C WE  WILL  USE  THE  APRAYS  ERRSAV.  TFS T 1 » AND  TEST2  IN  THE  TEST 

C FOR  COPRYCT  ORDER.  THE  SOLUTION  AT  LEVEL  K = 0 . 1 » 2 . 3 , 4 , OR 

C 5 SHOULD  BEHAVE  AT  ORDER  4 . 8 ♦ 1 2 * 1 6 . 20 ♦ OR  24  RESP.  WHEN 

C THE  MESH-SPACING  IS  HALVED  THE  ERROR  ESTIMATE  SHOULD  DECREASE 

C BY  FACTORS  OF  1 6 . 256 . 4096 . E T CETERA.  THE  ELEMENTS  OF 

C TEST  1 AND  TEST2  CONSTITUTE  (RATHER  COARSE)  LOWER  AND 

C UPPER  BOUNDS.  RESPECTIVELY.  ON  WHAT  IS  AN  ACCEPTABLE  LEVEL 

C OF  DECREASE.  WHILE  ERRSAV  SAVES  THE  ERRORS  FROM  THE  PREVIOUS 

c n to  allow  computation  of  these  factors  and  also  to  help 

C DETERMINE  the  LEVEL  AT  WHICH  INTERPOLATED  VALUES  ARE 

C ACCFPTED  when  N is  DOUHl  ED. 

DIMENSION  ERRSAV (6).TEST1(6).TEST2(6) 

DATA  TESTI/8. .64. »S12. .4096. .3.E4.2.E5/ 

DATA  TFST2/32. . 1024. ,3.E4, 1 ,E6. 3.E7, 1 .F9/ 

C WE  WILL  USE  LOCAL  LOGICAL  VARIABLES: 

C DIVNFW  (TRUE  WHEN  NEwTON^S  METHOD  HAS  DIVERGED) 

C NEWN  (TRUF  WHEN  WE  are  COMPUTING  AT  A NEW  VALUE  OF  N) 

C RKRROP  (TRUE  WHEN  WE  HAVE  DETECTED  DOMINANT  ROUNDING  ERROR) 

C CONVOH  (TRUE  WHEN  WE  have  HAD  CONVERGENCE  OF  NEWTONXS  METHOD 

C AT  LEAsT  ONE  TIME.) 

C NOEVM  (TRUE  WHEN  E AND  OE Y HAVE  ALREADY  BEEN  EVALUATED.) 

LOGICAL  Y I NT  .RERROR.DIVNEW.CONVOB.NEWN.NOEVAL 
COMMON  /flag/  IEI.AG 
COMMON  /K DUNT / K OUN T 1 • K DUN T 2 
DATA  Nw  A x /2R6/ 

DATA  N**  A v m / 1 28/ 

C CHECK  THF  VALIDITY  OF  ThT  INPUT  PARAMFTERS 


onooon  ooooo  n oonoori 
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IF(A  ,E0.  B)ST0P?6 

C FOR  THE  PRESENT.  WE  STOP  WHEN  THE  LINEAR  INTERPOLATION 

C PROBLEM  HAS  NO  SOLUTION. 

IE ( ( <ALEO*A*RV0) *ALEN)  .EQ.  ( ( AlEN°B.RVN) °ALE0) ) ST0P?7 
C 9.B9E-13  EQUALS  l.E-12  - EPS  FOR  OUR  PURPOSFS 

IEIT0L  .LE.  9.99E- 13)STOP30 
I E ( .NOT.  Y I N T ) N = 8 

INITIALIZE  SOME  VARIABLES 

THE  INITIAL  VALUE  OE  N IS  STORED  aS  NO  5 WHEN  N = N0  NO  BACKGROUND 
INFORMATION  EXISTS. 

MAX*  IS  THE  HIGHEST  PREVIOUS  LEVEL  OE  DEFERRED  CORRECTION  THAT 
HAS  BEEN  REACHED.  WHILE  K IS  THE  PRESENT  LEVEL. 

EPNEWT  1=  The  TOLERANCE  FOR  CONVERGENCE  IN  NEWTON^S  METHOD. 

? M A =0 
N0  = N 

E PSNw  T = TOI  /?  . 

RERROR=  .EalsE  . 

CONVOB=.  FALSE. 

K =0 

BEGINNING  OF  COMPUTATION  with  A new  value  OE  N. 

3 ERRoi  ' = 1 . E?0 
NP1 = N ♦ 1 
NP?=N.? 

NP3rNJ.3 

i-MAX  IS  MAXIMUM  ALLOWABLE  NUMBER  OF  EXTRAPOLATIONS  FOR  A 
GIVEN  N.  THERE  ARE  THREE  LIMITS.  THERE  MUST  BE  ENOUGH  POINTS 
TO  COMPUTE  THE  EXPANSIONS?  6 IS  DETERMINED  BY  THE  SIZE  OE  THE 
RELEVANT  ARRAYS;  M A XK ♦ ? IS  A RESTRICTION  ON  THE  NUMBER  OF 
EXTRAPOLATIONS  WITHOUT  HISTORY. 

KM AX=  < NP 1 -A ) /A 
KMAX  = MJN0 (KMAX»6*MAXK  *2) 

NE WNr . TRUE . 

E0 ( 1 i . ALFO 
E 1 ( 1 ) hPVO 
EO ( NP  3 ) =ALFN 
E 1 ( NP  3 ) =8 VN 
H — ( B - r ) /N 

WE  START  PREPARING  TO  USE  NE WTON-*S  METHOD  TO  SOI  VE  THE 
NON-i INEAR  EQUATIONS  APPROPRIATE  TO  THE  SPLINE  COLLOCATION  METhOO 
USED.  The  VARIOUS  CORRECTION  TERMS  MUST  BE  INITIALIZED 
APPROPRIATELY.  TO  ZERO  AT  LEVEL  ZERO. 

ITLIm  ANO  ITLIM?  LIMIT  THE  NUMBER  OE  NEWTON  ITERATIONS  AS 
DESCRIBED  LATER. 

I TL  I M t l 0 
I TL I M2  = 1 3 

» ( 1)  r ' 

» (NP  l i 

">0  10  I=?,N 
10  X ( I ) =X ( 1 ) ♦ ( I - 1 » 

DO  12  1=1. NP3 
YPC ( T ) =0 . 

YPCOt.O  ( I ) =0  . 

12  S( I ) =0. 

I E ( N .FQ.  NO ) GOTO  1 R 
C FILL  IN  Y VALUES  AETFR  DOUBLING  N 

CALL  GETVAL (N.H.RHSTOT . Y.YP.AYPP) 

C WE  USE  DIEEERENT  YP  VALUES  AFTER  AN  INTERPOLATION 
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C TO  EVALUATE.  1 HE  FUNCTION.  SO  EVALUATION  IS  REPEATED. 

OO  IS  1=1. NP1 

FU  ( I ) -zF  ( X ( I ) . Y ( I ) . YP  ( I ) ) 

e o ( i ♦ i ) = -or y ex ( i ) « y < i > * yp m > 
is  El (I ♦! ) =-OFYP ( X ( I) . Y < I ) , YP ( I ) ) 
xOUNT 1 =KOUNT 1 »NP 1 
K OUN  T ? = KOUNT  <?  »NP  1 
IF(k  ,FQ.  0 ) GOT  0 38 
CALL  M A THAN < W.RHS, EO, E 1 . WA .N.H. 1 ) 

IF  ( I F i AG  .GT.  0 ) GO  TO  100? 

CALL  PRSUR (N.?8 . ERRNOR . X . RHST OT ) 

100?  CONTINUE 

C NFXT,  EVALUATE  The  NEW  CORRECTIONS  S AND  YPC  FOR  THE  NON-LINEAR 

EQUATION^  that  MUST  re  SOLVED  AT  The  new  level  k after  an 
INTERPOLATION  FOR  Tor  VALUES  ACCEPTED  AT  LEVEL  K-l. 

WHEN  NEWN  IS  TRUE  . OF Y IS  EVAi  UATED  DURING  THE  NEWTON 
ITERATION . 

NEWN=. FALSE. 

ITL I M=  T 
I Tl IM  j = 5 

CALL  SPLDCGU  .N.FU.S,  IFRROR) 

CALL  0 YDC  G ( K . N . H » F . Y PC ♦ I ERROR ) 

DO  17  1 = 1. NP1 

1 7 YPCOl  0(1) =YPC ( I ) 

PO  is  1 = 1. K 

18  ERRSAV  ( I ) =FRRSAV  ( 1 ) • 1 G . <>  <M  - I ) 

GOTO  JR 

19  IF  ( V ; ;T  ) GOTO  3? 

OO  ?0  1=1. NP3 

?0  R^ST  'T  ( I ) =0 .0 

IF  < ( » VO  ,NF.  0.)  .OR.  («VN  .NF.  0 . ) 1 GOTO  38 
CALL  LINST (GO.ALFO.GN.ALFN.A.B.N.H.RHSTOT ) 

GOTO  Fq 

3?  CALL  GFTST ( W.RHSTOT ,WA .N.H) 

IFdFLAG  .EQ.  0 ) C ALL  PR  SUB  ( N . ?h  . T WOP  I . X . RhS  TO  T ) 

C m a I N PORTION  OF  NEWTON^S  METHOD 

C WE  ENTER  HERE  WITH  an  INITIAL  GUESS  FOR  Y.  EITHER  THE  LINEAR 

C INTERPOLATION  OF  THE  BOUNDARY  DATA  OR  AN  IMPROVED  ESTIMATE 

C FROM  A LOWER  LFVEL  AT  THIS  N OR  A HIGHER  LEVEL  AT  A SMALLFR  N 

C OR  A USER  SUPPLIED  INITIAL  GUFSS. 

C AT  LEVEL  *=0  WE  ALLOW  11  ITERATIONS  BEFORE  TERMINATING  T HF 

C PROCEDURE.  AT  HIGHER  LEVELS.  ONLY  u ITERATIONS  ARE  ALLOWEn. 

C AT  LEVEL  K = 0 NEWN  IS  TRUE  AND  WE  EVALUATE  DFY  AT  EACH 

C ITERATE.  FOLLOWING  AN  INTERPOLATION  ACCEPTED  AT  HIGH  ORDER  F . DF Y . 

C AND  DFYR  ARE  EVALUATED  BEFORE  ENTERING  THE  LOOP  AND  ARE  NOT 

C RF  E V At  DATED  ON  ENTRY.  OF  Y AND  DE YP  ARE  NEVER  EVALUATED 

C WITHIN  THE  LOOP  IN  This  CASE. 

C NORMAi.  EXIT  OF  THE  NEWTON  ITERATION  OCCURS  WHEN  THE  NORM 

C (RNORM)  OF  THE  CHANGE  IN  THE  SOLUTION  Y IS  LESS  THAN  OR 

C EQUAL  TO  HALF  THE  USFRhS  TOLERANCF  TIMES  THE  NORM  OE  Y. 

C ABNORMAL  EXITS: 

C (1)  THE  NORM  OF  THf  RFSIOUAi  INCREASES  FROM  STEP  N TO 

C N ♦ 1.  WHERE  N > 0. 

C OR 

C ( ? ) THE  NIJMRER  OF  ITERATIONS  EXCEEDS  ITlIM  AND  CONVERGENCE 

C IS  NOT  EXPFCTED  RY  ITLIm^. 

IF  THF  CHANGE  IN  Y (RNORM)  INCREASES  BUT  THE  RESIDUAL  DECREASES, 
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THIS  IS  TA-sn  AS  A SIGN  THAT  ROUNDING  ERROR  IS  PREVENTING 
CONVERGENCE  TO  THE  TOLERANCE  AND  THE  PROGRAM  GIVES  AN  ERROR 
MESSAGE  AND  TERMINATES. 
lMr>EX=0 

DIVNEW=. FALSE. 

NOEVAL= . TRUE . 

I E ( ( N .EQ.  NO)  .OR.  (x  .GT.  0 ) ) NOE V AL = . E ALSE . 

CALL  GETVAL (N.H.RHSTOT.Y,YP,AY?P) 

CALL  BC ( ALE  0. BV0 t GO »ALEN»BVN*GN»H»N»RHSt RHSTOT) 

RE  SN  = AM  A X 1 ( ABS ( RHS ( 1 ) ) » APS  i RhS (NP  3 ) ) ) 

IE (NOEVAL ) GOTO  6S 
O 60  1=1. UP1 
YPTEMP=YP ( I ) -YPC ( I ) 

EU ( I ) =E ( X ( I ) »Y ( I ) « Y P T E M P ) 

IE ( .NOT.  NEWN ) GOTO  60 

E0 ( I ♦ 1 ) =-DEY  (X ( I ) .Y ( I ) .YPTEMP) 

E 1 ( I ♦ 1 ) =-DEYP ( X ( I ) . Y ( I ) , YPTEMP) 

CONI INUE 

KOUNV 1 =K OUN  T 1 *NP1 
I E ( .NOT.  NEWNIGOTO  65 
K0UNT2=K0UNT2*NP1 
NOEVAL  = .EAL  SE. 

DO  70  1=1. NP1 

RHS ( I ♦ 1 ) =S ( I ♦ 1) *EU ( I > -AY2P  < I > 

IE (RE  SN  .LT.  ABSIRHS ( I ♦ 2 ) ) )RESN  = ABS(RHS ( I ♦ 1 > ) 

CONTINUE 

r iow  = P 

IE (NE-N) I JOB=0 

CALL  MATMANt  V.RHS.E0.E1  .WK.N.H.  I JOB) 

RNORh=0 , 

RTNORM=0. 

HO  100  1 = 1.  »P3 

RHS  ToT ( I ) =RHST0T ( I ) *RHS ( I ) 

I E ( RNORM  ,LT.  ABS (RHS ( I ) ) ) RNORM  = A0S (RHS ( I ) ) 

I E ( R TNORM  .LT.  ABS  (RHSTOT  ( I ) ) ) RTN0RM.-A8S  (RHSTOT  ( I ) ) 

AYNORM  IS  AN  APPROXIMATION  TO  THE  NORM  OE  Y BASED  ON  THE 
THEORETICAL  BOUND  ON  THE  RATIO  OE  THE  NORMS  OE  THE  SPLINE  COEE. 
AND  THE  FUNCTION  Y. 

AYNORM=RTNORM/3. 

INDEX= INDEX ♦ 1 

IE ( INDEX  .EO.  1 )RNORMl=RNORM 
EPS=EPSNWT "A yNORM 

IEIIELAG  .EO.  0)PRINT  470  , INDEX. RESN. RNORM 

I E ( RNORM  .LE.  EPSIGOTO  3B0 

IE ( INDEX  .LE.  2)  GOTO  120 

IEfRESN  .GT.  RSNSAVIGOTO  300 

IF ((RNORM  .GT.  RNc  A V ) .AND.  CONVOB>GOTO  350 

RNS2--RNS  A V 

RNS  AV  = RNORM 

RSNS  A V = RE  SN 

IE ( INDEX  .LE.  ITLIMIGOTO  40 
IEIINDEX  .GT.  I T L IM2 ) GOTO  300 

I E ( ( RNORM  » (RN0Rm/RNS2)*« ( ITL IM2-INDEX) ) .LT.  EPS)GOTO  40 
NFWTON-.S  METHOD  IS  DIVERGING  AND  WE  MUST  DOUBLE  N IE 
POSSIBLE.  IE  THE  EINAL  NEWTON  ITERATE  SOLVES  THE  NON- 
I INEAR  SYSTEM  BF  T TER  THAN  THE  FIRST  ITERATE.  THEN  WE  GET 
an  ERROR  ESTIMATE  FOR  USE  IN  TESTING  ORDER  AND  WE  USE 
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V At  i if  s 1 N T E RP0L  A TE  D FROM  ThE  LAST  ITERATE  IN  ORDER  TO 

START  NEWT  0N*S  method  at  level  zero  and  with  n doubled  I 

IE  THE  FINAL  ITERATE  WAS  WOPSE * HOWEVER*  WE  START  WITH 
THE  LINEAR  INTERPOLATION  OF  THE  BOUNDARY  DATA, 
n I VNE W= .TRUE  . 

I F ( ( RNORM 1 .OT.  RnORM ) .OR.  CON V OB) GOTO  400 
IF / N ,GT.  NMAXHIGOTO  p00 
IE (YlNT)GOTO  R?0 
N = J * N 
GOTO  ? 

WE  NOW  THlNrv  THAT  ROUNDING  ERROR  HAS  CONTAMINATED  OUR 
RESULTS,  SO  WE  PREPARE  TO  OUIT  COMPLETELY  AFTER  ESTIMATING 

THE  ERROR  we  DID  ACHIEVE. 

REPROR= . TRUF . 

GO t 0 AO  0 

WE  NOTE  FOR  FUTURE  REFERENCE  THAT  CONVERGENCE  HAS  BEEN 
ORTAINEO  AT  LEAST  ONCE. 

CO-  VOB=.TRUE. 

we  WANT  TO  ESTIMATE  THE  DISCRETIZATION  ERROR*  ESSENTIALLY 
BY  ONE  STEP  OF  NEWTON-S  METHOD.  WE  SET  UP  THE  RIGHT  HAND 
SIDE  RHS  NEEDED  FOR  THE  LINEAR  EQUATIONS  OF  ThIS  ONE  STEP « 

WH 1 1 F COMPUTING  RHS  wF  WILL  ALSO  COMPUTE  S AND  YP,  THE 
CORRECTION  (ESTIMATING  LOCAL  TRUNCATION  ERROR)  THAT  WE  WILL 
NEED  FOR  THE  NON-LINFAR  EQUATIONS  AT  THE  NEXT  LEVEL  UP  IE 
WE  Gt T THAT  FAR. 

A ♦ 1 
IT,  I M=  3 
IT  I M?=5 

CALL  SPLDCG(K,N.EU*RHS, IEHROR) 

Call  DYDCGUa  ♦N.h,EU*YPC*  I ERROR  ) 

CALL  BC(ALFO*BVO*GO,ALEN*BVN,GN*H,N,RHS*RHSTOT) 

^0  450  I -i.  * NP2 
ST  =PHS(I) 

RHS ( I ) =S  < I > -STF 
S r T ) =STF 

RHS ( I ) =RHS ( I ) -E I ( I ) * ( YPC ( 1-1 ) - YPCOL  0(1-1  ) ) 

YPCOLDi 1-1 ) = YPC ( 1-1  ) 

WITH  THE  NEW  RIGHT  HAND  SlDF*  WE  GENERATE  AND  SOLVE  THE 
NEW  SYSTEM. 

CALL  MATMAN(W*RhS*E0*E1*WK,n,H,?) 

ERROR  CONTROL  AND  DECISION  CENTER 

Rhc,  NOW  CONTAINS  OUR  ESTIMATE  OE  DISCRETIZATION  ERROR. 

WE  COMPUTE  ITS  NORM  IN  ERRNOR.  ALLOWING  US  TO  TEST  E OR 
CONVERGENCE  OE  THE  OVERALL  ALGORITHM. 

NEWN=. FALSE  . 
t «p-.iOR  = 0. 

Y NORM  r 0 . 

GO  460  1=1. NP1 

T E = ABS  ( (RHS  ( I > ♦ 4 . R H S (1*1)  .RHS  ( I *2)  >/6.  ) 

I E ( TE  .GT.  ERRNOR) EPRNOR=TE 

IFtYNOPM  .|,T.  ABS  ( Y ( I)  ) ) YNORM  = ABS  (Y  ( I ) ) 

IE ( I E L A G .GT.  0 ) GOTO  1001 
KUk-1 

print  4Z0f INDE * .RESN.PNORM 
CALL  PRSUtHN.f  1 .ERRNOR  . X .RmSTOT  ) 

CONTINUE 
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C NEXT  WF  TEST  OUR  ESTIMATE  OF  OVERALL  ERROR  AGAINST 

The  USER^S  TOLERANCE.  tPPNOR  ESTIMATES  The 
DISCRETIZATION  EPPOR  FOR  ThF  EXACT  SOLUTION  OF  THE 
NON-LINEAH  SYSTFM,  while  PNOPM  ESTIMATES  The  ERROR  IN 
SOLVING  THAT  NON-LINFAP  SYSTEM.  ERREST  ESTIMATES  THE 
TOTAL  ERROR  UNDER  THE  PESSIMISTIC  ASSUMPTION  THAT 
THE  OTHER  ERRORS  SUM. 

E PWEST=ERPNOR*RNORm 
T 01  LOC=TOL*YNORM 
TFST  FUR  SUCCESS 
ifieprest  .LT.  TOLLOOGOTO  QOO 
TERMINATE  algorithm  if  rerror  true 
IF (PFRROR)GOTO  080 

if  k ,lf.  maxk  (so  that  there  is  history  to  compare 

AGAINST)  WE  NEXT  CHECK  WHETHER  THE  ERROR  IS  BEHAVING 
AT  THE  CORRECT  ORDER. 

1F(K  .GT.  MAXKJGOTO  480 
RAT=ERRSAV(K) /eppnop 
ERRSAV (K ) =ERRN0R 
IF ( I FL  AG  .EQ.  0 ) PR  I N T 4 78  » R A T 
THIS  IS  THE  TEST  FOR  CORRECT  ORDER. 

IE (PAT  .LT.  TEST  I ( K ) .OR.  PAT  .GT.  TEST2 ( K ) ) GOTO  700 
480  EPRSAV (K ) =EPRNOP 

THF  NEXT  STATEMENT  TESTS  WHETHER  THE  PREVIOUS  DEFERRED 
CORRECTION  IMPROVED  THE  ACCURACY  ACCEPT  ABl  Y AND 
WHETHER  AN  ERROR  ESTIMATE  WOULD  BE  ALLOWED  AFTER 
THE  NEXT  CORRECTION. 

IFiEPRNOR  .GE.  .l*ERPOLD  .OP.  K * 1 .GT.  KMAX ) GOTO  700 
CHECK  FOR  NEW1 ON  D I VEPGENOE 
IF ( O I VN£W ) GOTO  700 
E RROl n-EPRNOR 

GO  AHEAD  AND  SO!  VE  THE  SYSTEM  AT  THE  HIGHER  ORDER 

GO T 0 08 

WE  NFFD  TO  INCREASE  N AND  THE  NEXT  SECTION  OF  CODE  HANDIES 
THAT.  WF,  FIRST  COMPUTE  THE  LEVEL  INTERPOLATED  VALUES  WILL  BE 
ACCEPTED  AT.  THEN  WE  DO  THE  INTERPOLATION. 

700  I F ( N .GT.  NMAXH)GOTO  800 

M A X K =K 
NOL  n=N 
N = ?*N 
K I NT  rK 

ERRO!  0 = AM  I N 1 (FRROLD.ERRNOR) 

C DECIDE  LEVEL  TO  INTERPOLATE  TO 

DO  710  I — 1 , K 

I F ( EPPOLD  .LT.  ERRSAV  ( I)  *>  ( -I  )>  GOTO  710 

C-DTO  /?0 
710  CONTIMUF 

720  K = I - 1 

I F < o i NE  W ) K = 0 

C INSPI  PERFORMS  THE  INTERPOLATION  AND  FILLS  IN  Y VALUES. 

CALL  INSPL (KINT.N«Y*w»S»RhSTOT»WK»K.H> 

GOTO  1 

C THIS  IS  ThE  N TOO  BIG  EXIT 

800  PRINT  8 1 S 
GOTO  900 

C 820  IS  ThE  ERROR  EXIT  TAnEN  WHEN  AN  INITIAL  GUESS  FAILS. 

820  PRINT  8?s 
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CALL  SPLGFN (w.N.H.RHSTOT ) 

RETURN 

C 980  is  thF  R0UN0  I MG  ERROR  FXIT. 

9*0  PRIM  QQO 
r, QTO  Ron 

4 70  FORMAT U l 0 « 2f  20 . 1 0 > 

875  FORMAT (ROX .F 12.4) 

915  FORMAT  i * M TOO  BIG*) 

825  FORMAT!*  ILLEGAL  TO  REINITIALIZE  WHEN  INITIAL  GUESS  GIVEN*) 
920  FORMAT!*  YNORm  = *,E20.10> 

9R0  FORMAT!*  ROUNDING  ERROR  EXIT*) 

ENO 


SUBROUTINE  GFTVAI  (N.H.RhS.Y, YP, AY2P) 

this  ROUTINE  Takes  the  valUF  OF  THE  SPLINE  COEFFICIENTS 
AND  COMPUTES  THE  VALUES  OF  Y.YP.  AND  AY2P  FOR  USE 
IN  EVALUATING  The  USER  FUNCTIONS  F.DFY.AND  DFYP. 

DIMENSION  RHS  ( 1 ) . Y ( 1 ) . YP  { 1 ) , AY2P ( 1 ) 

NP I -N  ♦ 1 
00  30  Irl.NPl 

Y ( I ) = (PHS ( I ) ♦A.ORHS ( I ♦ 1 ) *RHS ( I *2) ) /6, 

30  YP  U ) = (RHS ( I *2) -RHS( I ) > / (H*2. ) 

AY2P(1)=7.*RHS(1)-16.5*RhS(?) ♦ 14.*RHS (3)-7.*RHS(4) 

1 ♦ 3 . *RHS (5) -,5*RHS (6) 

DO  40  1=2, N 

AY2P ( I ) = . 5*RHS ( 1-1 ) ♦4.*RHS< I ) -9.*RhS< 1*1 ) ♦4,*RHS( 1*2) 

1 ♦ .5«R-iS  < I ♦ 3) 

40  CONTINUE 

AY  2P ( NP 1 ) =-.5*PHS (N-2) *3.*RHS (N-l ) -7.*RHS (N)  ♦ 14.*RHS (NP1 > 
1 -16.5*PhS(N*2) ♦7.*RhS(N*3) 

DO  50  1=1. MP 1 

50  AY2P ( I ) =AY2P ( I ) / <6.*H*H) 

RETURN 

END 


SUBROUT  TNF  M A TM AN ( W , RHS , EO ♦ F 1 « WK , N , H * I JOB) 

THIS  ROUTINE  GENERATES  THE  SOLUTION  MATRIX  W FROM  THE 
VECTORS  OF  PARTIAL  OFRIVATIVE  VALUES  EO  AND  El. 

OIV1  AND  0IV2  ARE  ASSUMED  P«ESERvFD  BETWEEN  CALLS  TO 
FACILITATE  REPEATED  CALLS  WITH  ThF  SAME  EO  ANO  El 
VALUES  EQR  DIFFERENT  RIGHT  HAND  SIDES.  I JOB  IS 
DEFINED  IS  IN  LEOT1R. 

AFTER  Thf  MATRIX  IS  GENERATED  The  LINEAR  SYSTEM 
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C IS  SOL  VF  0 BY  CALLING  LEQT1R. 

DIMENSION  W (250,5) ,RhS  (250)  ,E  0 ( 250)  ,E  1 ( 359)  ,WK  ( 777) 

HSQ=H*H 
NP 1 =N ♦ 1 
NP2  = N ♦ ? 

NP3  = N ♦ 3 

HSQ6  I - 1 . / ( 6 . °HltH  ) 

IFUJOH  .FQ.  2 ) GO  T 0 100 

C POWS  1 AND  N ♦ 3 DFPEND  ONLY  ON  THE  BOUNDARY  DATA  WHICH  WE  STORE 

C IN  EO  AND  El. 

W(1.3)=(F0<1) »H-3. *E 1 ( 1 ) ) / ( 6."H> 
w ( 1 ,4)  =2.<>E0  ( ! ) /3. 

W(1.5>  = (rO(l)^H*3.<»El(l))/(6.<>H) 

w ( NP  3 * 1 )=(E0 (NP3)*H-3.ftE 1 (NP3) ) / ( 6 . *H ) 
w (NP3. 2) -2.  »F0 ( NP 3 ) / 3 . 

w (NP3.3)  = (EO (NP3) °H*3.*F1 ( NP  3 ) >/<6.*H> 

C ROWS  3 THROUGH  N-l  ARE  STANDARD  FIVE  BAND  ROWS  BASED  ON  THE 

C DIFFERENTIAL  EQUATION. 

DO  50  I = 3 » NP 1 

CON  1-3 . °H°E 1 ( I ) 

COnB-hcq*fO  ( ! ) 

W(I,2)=(4.-C0N1*C0N2> *HSQ6 I 
W <1 , 3)  = (-0. *4. »CON2) *HSQ6I 
W ( I .4) = <4. *C0N1 *C0N2) «HSQbI 
W <1  » 1 ) =0.5*HSO6I 
50  W ( I » 5 ) =0 . 5*hS06 I 

C ROWS  2 AND  N*2  HAVE  SIX  ELFMENTS  IN  THE  ORIGINAL  MATRIX;  PRE- 

C ELIMINATION  IS  USED  TO  ELIMINATE  TWO  ELEMENTS  OF  EACH. 

CON  1 = 3 . »E  1 ( P ) -“-H 
C0N2=HSQ»E0  (?) 

DI V2=-2.0* < 3. *W  <4 ,4) /HSQ6I ) 

W (2.2) = ( 7. -CONI -C0N2) «HSQ6I *DI V2*W (3*1) 

W (2. 3)  = (- 1 6.5*4.*C0N2) °HSQ6I *DI V2»W (3.2) ♦ W < 4 , 1 ) 

W (2 .4) = ( 14. ♦ CONI *C0N2>  *HSQ6I *DI V2»W <3,3)*W(4,2> 

W (2»5) =-7.*HS06I *DI V?*W (3*4) *W (4,3) 

CON  1 =3  . »F 1 ( NP2 ) * H 
CON2  = HSQi>EO  ( NP2 ) 

D I VN  = -2 . A * ( 3. *W (N*2) /HSQ6I ) 

W (NP2* 1 ) = - 7 , «HSQ6 1 *DIVN*W(NP1,2)*W(N,3) 

W (NP?*2) = ( 14. -CONI *CON2 ) *HSQ6 I ♦ D I VN* W ( NP 1 . 3) *W (N*4> 

W ( NP?  * 3 ) = (- 16.5*4. *C0N2) *HSQ6I *DlVN*W (NP1 *4) *W (N*5) 

W (NP2  *4 ) = ( 7 . *CON 1 *CON2 ) *HSQ6I *DI VN°W (NP1 *5) 

100  CONTINUE 

IF ( I JOB  ,FO.  l)GOTO  150 

PHS (2) =RHS (2) *RHS(4) *DIV?«RhS(3) 

RHS  ( NP2  ) =RHS  (NP2)  ♦ RHS  ( N ) *DIVN<>RHS  (NP1  ) 

150  CALL  LFOT1B(W*NP3*2*2*Rhs,wk, I JOB, IER) 

IF ( IFR  ,nF • 0 ) S T 0P44 

RF  TURN 

END 
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SUHPOUT INF  INSPL  (K,N*Y,W»S»RHS*WK,KNFW*H) 

C THE  difference  correction  GENERATION  part  of  the  next 

C ThREF  ROUTINFS  IS  BASED  ON  U2DCG  OF  PFREYRA. 

C SEE  r H 1 0 H ORDER  FINITE  DIFFERENCE  SOLUTION  OF 

C DIFFERENTIA!  EUUATIONSE.  BY  V.  PFREYRA,  STANFORD 

C UNIV.  CS  REPORT  STAN-CS-73-348»APRIL , 1973. 

C THIS  ROUTINE  DOES  THE  INTERPOLATION  OF  THE  Y VALUES 

C AND  THEN  CALLS  GETCO  TO  GET  SPLINE  COEFFICIENTS. 

DIMFNMON  Y(1),W(1)»S(1>  , RHS  ( 1 ) »WK  (1) 

DIMENSION  A (28) ,C (28)  ,CCON (28)  , AEXP  (28) 

DATA  CCOn/2*0. .2. .0. * -2. . 0. . 2. * 0., 11. 3333333333333*0. »-l 98. ,0. 

1 2157.  ,0.,  -12131. 3133333333,0  . .-3.4  82286666666666E5 ,0.  , 

2 19492202. *0. .-5. 9*299 148666666E 8,0 . , 7 . 762694  1 35  33333E  9 , 0 . , 

3 6.07798554731 133E1 1 . 0 . , -6 . 8422 1 4 7682 7 757E 1 3 . 0 . / 

DATA  A/1 .0*23*0.0/ 

N02=N/2 
N02 1 =NO? ♦ 1 
DO  10  1=1.5021 
10  W ( I ) = Y ( I ) 

KK 1 = 4«K 
I I =N02 1 -K K 1 

KK  =KK 1-1 

KM  I D= ( KK 1 ♦ 1 )/? 

KMID1 =KM ID” 1 
KK 1 P 1 =KK 1 ♦ 1 
KK I JP1 =KK ♦! I ♦ 1 

c asymmetric  approximations  near  left  boundary 

DO  20  1=1 .KM iol 
XI=I*0.5 

CALL  COEGEN (KK 1 , x I ,C, A ) 

C Uc  E RFFiECTION  NEAR  RIGHT  BOUNDARY 

AC?'  ->  = 0 . 

DO  12  J = 1 « K K 1 

12  ACIJM  = ACUM*C  (KK1P1-J)  *Y  ( I I ♦ J) 

IF(K.NE.O)  S ( KK I I P 1 - I ) = ACUM 
IF(K.EO.O)  S (N021 • 1 ) =ACUM 

AC1 1^  = 0 

DO  15  J = 1 , K K 1 

15  A CUM = A CUM ♦ C ( J ) * V ( J ) 

20  S ( I ) = ACUM 

C SYMMETRIC  APPROXIMATIONS  CONSTANT  THROUGHOUT  CENTER  RANGE 

X I=KMID*0.5 

CALL  COEGEN(KK 1 , x I ,C, A ) 

NF  = N02 1 -KM  I D 
DO  40  I = ► MlD.NF 
AC"  ‘ = 0. 

I I = I -KM  I n 
DO  38  J= 1 « KK 1 

38  ACUM  = ACl)M*C  ( J)  <-Y  ( I I ♦ J) 

40  S ( T ) = ACUM 

DO  60  1 = 1,  '102 
12  = 2*  I 

Y ( 12-1  ) r ( I ) 

60  Y ( 12) =S ( n 

Y (N* 1 ) = w (N021 ) 

NP 1 =N  ♦ 1 
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DO  7 0 1 = 1.  NP 1 
70  -7hS  < I * 1 ) = Y ( I > 

Call  GETCO  < * ,RHS* WK ,N,H,KNEW ) 

RE  T URN 

END 


SUBROUT  INF  SPLDCGtK.N, Y,S, I ERROR) 

T H I s ROUTINE  COMPUTES  THE  DIFFERENCE  CORRECTION  FOR 
THE  BASIC  METHOD.  IT  COMPUTES  THE  CORRECTION  FOR 
POINTS  l....,N*l.  SEE  THE  REFERENCES  ON  DIFFERENCE 
CORRECTION  FOR  this  method. 

DIMENSION  A (20) , Y < 1 ) ,S  < 1 ) ,C (28) ,GA (28) ,GB (28) 

DATA  A/4<>0  . , 6.66666666666667F-2.0.  , 1 . I 904  76 1 904  76  I 9E  - 1 » 0 . , 

1 -1 .08888888888b88,0.*8.  1 8 1 8 1 8 1 8 1 8 1 8 1 8 * 0 . , - 3 3 . 3223443223443 , 0 . , 

2 -725.472222222222,0. ,-3185.00065359477,0. ,-7.8 1 97256 140 3509F5, 

3 0 . , 8 . 40 1 1 84 131 31 31 3E6»0.» 5.505421 69 1 425 1 2E 8 , 0 . , 

4 -5. 26 3242 1 29444 14E 10,0. ♦ 2 . 66683028642333E 1 2 , 0 . / 

DATA  GA/40O., -1.9333333333333 3, -10. , - 34 . 880952 3809524 ,- 1 05 . , 

1  -285.755555555555,-714. , - 1 . 75 1 8 1 8 1 8 1 8 1 8 1 8E 3 , -4*20 . , 

? -i . 3 123322 344 3223E 4, -3 7466. ,-9. 7746638888888E4 , 186777.5, 

3 -2. 147066601 3071 8E5,-1 . 544223333333 33E6 ,- 1 . 4 73972056 1 40 35E 7 , 

4 -51816952. , 7. 52986 0079 79798F 7, 1 748890 1 50 ., 6 . 90 1 3763 1 0809 1 7E9 , 
5-5.1 1509895563333F 10, -6.861762420 12174E 11 *-5.0 093686 1343333F 11 , 
6 34653665345130. , 1 999704250626220 . / 

DATA  GB/4»0. .-1 .93333333333333, 1 0 .,- 34 . 0809523809524 , 105. , 

1 -285.755555555555,714. .-1 . 75 1 8 1 8 1 8 1 8 1 8 1 8E 3 , 4620 . , 

2 - 1 . 31 23322344322 3F4.37466. , -9. 7746638888888F4,- 1 867  77.5, 

3 -2. 14706660 130 7 18E5, 1 . 544223 3 3 333333F 6 , - 1 . 4 7 3972056 1 40 35F 7 , 

4 51816952. ,7.52906007979798^7,-1748890150. ,6. 90 1376 3 10809 17F9, 

5 5.1 1509M9556  3 333F 1 0 , -6 . 86 1 762420 1 2 1 74f  1 1 , 5 . 0 0 93606 1 34  3333E 1 1 , 

6 346536653451 30. ,-1999784250626220./ 

C CHEC*  INPUTS 

I F ( ( k 0)  .OR.  (K  .GT.  (N-3)  /4)  ) GOTO  100 

KK 1 =4* (K  ♦ 1 ) 

KK=KK 1 -1 

KMI()=  (KK  1 ♦ 1 ) /2 
If  RROR  = 0 

KM  I ' 1 =KM ID-  1 
K l NT  tKK  1 
I I =N-KK 
KK 1P1 =KK 1 ♦ 1 
KK I IR3  = K<  , r I , 3 

C SPECIAL  APPROXIMATION  AT  LEFT  BOUNDARY 

CALL  COEGEN (KK 1 , 1 . ,C»GA ) 

ACM'  =0. 

DO  2 J=  1 . KK  1 
2 ACUM=ACUM*C ( J) »Y ( J) 

S ( 2) = ACU * 

C ASYMMETRIC  APPROXIMATIONS  NEAR  LEFT  BOUNDARY 

DO  5 I=2.KMJD1 
XI  = I 

call  COEGENfKKl ,xl ,C»A) 
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U';*-  REFLECTION  NEAP  RIGHT  BOUNDARY 
ACUM=0  . 

DO  3 J= 1 » KK  i 

ACUM  = ACUM*C  (KK  IP  1 - J ) «•  Y ( I I ♦ J ) 

S (KK I I P 3- I ) =ACUM 
ACM  4 = 0 

DO  4 J=  1 * KK  1 
ACUM=ACUM*C  < J ) •»>  Y ( J ) 

S ( I ♦ ] ) =ACUM 

symmetric  approximations  constant  throughout  center  range 

K I N T =KK 
X I =h.M  I D 

CALL  COEGFN (K INT «X  I tC» A ) 

NF  = N* i -K INT  *KM  I D 
DO  40  I=*mID.NF 
AC  m=0. 

I I=I-KMin 
DO  38  J= 1 • K I N T 
AC'IM  = ACUM*C  ( J ) <■•  Y ( I I ♦ J ) 

S ( t ♦ 1 ) = ACUM 
x I=KK  1 
I I =N-KK 

CALL  COFGFN (KKl , x I ,C.GB) 

ACUM  = 0 . 
nO  60  J=1 .KKl 
60  ACUM  = ACUM*C  < J)  <>Y  ( I I ♦ J) 

S ( N ♦ 2 ) =ACUm 
RFTURnI 
100  I ERROR  - l 
RT TURN 
END 


SUBROUTINE  DYDCGD-  »N»-m  Y.S*  IERROR) 

C THIS  ROUTINE  COMPUTES  THE  CORRECTION  TERM  FOR  YP  IN  THIS 

C METHOD. 

DIMENSION  DCON  ( a 7 ) .0<27)*C(27)*Y(1)*S(1) 

DATA  DCON/ 3*0. .-3. 3 3 3333333333 3 3F -2.0. .7. 93650 793650 794F -2,0. , 

1 -.  194444444444444,0.  . - 9 . 09  090  9 090909  09E -2 . 0 . . 1 2 . 645299 1 45299  1 . 

2 0..-226. 777777777778*0. .2677. 95506535948.0. .8104. 68810916179. 

3 0. . -2. 18828425555556^6, 0. . 1 . 02 324 1 1 75 1 2077F8. 0 . . 

4 -2.802S1085196407F9.0. . -2 . 7528 1 37 78 1 666 7F 1 0 . 0 . / 

C CHECK  INPUTS 

IF ( ( K ,l  ■ . 0)  .OR.  (K  .GT.  (N-3) /4) ) GOTO  100 
KK  1 =4«-  (K  ♦ 1 ) 

K K = K K 1-1 

KM  I 0= ( KK 1 ♦ 1 ) /2 
If  RROR  = 0 
KM  ID  1 :KM JO- 1 
I I =N ♦ 1 -KK 
KK 1 I I =KK1  ♦ I I 

no  2 i = i , kk 

2 DTI) =H*OCON ( I ) 
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AS-MM^TPIC  APPROXIMATIONS  AT  OR  NEAR  LEFT  BOUNDARY 
00  5 1-1  . k M 1 0 1 
XI  = I 

CALL  COEGEN(KK.XI.C«D> 

U-E  REFLECTION  NEAP  RIGHT  BOUNDARY 
ACUM  = 0 . 

DO  3 J-WKK 

3 ACUM=ACUM-C(KKl-J)*Y(Tr*J) 

S(KK1 I I-I  > = ACUM 

ACI  ."  = 0 
00  J=1,KK 

4 AC"M  = ACUM*C  ( J)  Y ( J) 

5 S ( I ) - ACUM 

SYMMETRIC  APPROXIMATIONS  CONSTANT  THROUGHOUT  CENTER  RANGE 
X I = n M I 0 

CALL  COEGFN(KK,xi,C.D) 

NE  =N ♦ 1-KK*KMI0 
00  90  1=  MIO.NF 
AC"^  = 0 . 

I I=I-KMIP 
00  3P  J=1.KK 

3S  ACUM=ACUM*C  ( J)  *Y  ( I I ♦ J) 

90  S ( I ) = ACUM 

RETURN 
100  IERR0R=1 

return 

END 


SUBROUT  t NF  SPl  GEN  ( W » N » H « RHS  ) 

THIS  ROUTINE  TAKES  THF  SPLINE  COEFFICIENTS  USED 
INTERNALLY  AND  COMPUTES  THE  RESULT  GIVEN  TO  THE  USER. 
DIMENSION  W (PS9.5) .RHS (259) 

N P 1 = N ♦ 1 
DO  30  ]-  1 • NP  1 

W ( I . 1 ) = (RHS ( I ) *9.oRHS (1*1) *RHS < I *2) ) /6. 

W < I . 2 ) = (RHS ( I *2) -RHS( I ) ) / (?.*H) 

W ( I . 3) = (RHS ( I ) -2 • *RHS (1*1) *RHS ( I *2) ) / (H*H> 

IF ( I ,E0.  NP1 ) GOTO  90 

30  W ( I .9)  = (RHS ( I *3) -RHS ( I ) *3.* (RHS  < 1 ♦ 1 ) -RHS ( I *2)  ) ) / <B.*H««3) 
90  RETURN 

END 


SUBROUTINE  COEGEN <N.XNP.C*RR) 

DIMENSION  C ( 1 ) *RR  < 1 ) 

oooo  OOOH  ooooo 

THIS  is  a SLIGHTLT  MOOIE IFO  FORTRAN  VERSION  OF  THE  ALGOL 
PROCEOtiRF  PVANO.  P.  90 1 OF 
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C SOLUTION  OF  vandfrmondf  SYSTEMS  OF  equations  BY 

C A.  BJOBCK  AND  V.  PEREYRA,  MATh.  COMP.  24.  PP. 891-903 

C WHERE  A DESCRIPTION  OF  THE  METHOD  CAN  BE  FOUNO. 

C 0-00  ooooo 

C IN  THIS  IMPLEMENTATION.  GIVEN  GRID  POINTS  1.....N* 

C SEPARATED  BY  A CONSTANT  H,  THE  COEFFICIENTS  ARE  COMPUTED  FOR 

C A POINT  (XNP-  ( XNP  J ) *H  BEYOND  GRID  POINT  [ LNP  ) . 

DO  1 I = 1 . N 

1 C ( I ) =RB< I ) 

NN=N- 1 
DO  6 1=1 , NN 
LL=N- I 
ALF I=I-XNP 
UO  h J- 1 ,LL 
K=N- J* I 

h C(K)=C(K)-ALFI»C(K-1) 

DO  8 1=1 ,NN 
K=N-  I 
X K I N= 1 • /K 
KP 1 =K ♦ 1 
DO  n J = K P 1 . N 

C(J)=C(J)*XKIN 

8 C U-l  ) =C  ( J-l  )-C  ( J) 

RETURN 

END 
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SUBROUTINE  -ETST (W.RHS.WK ,N.H) 

THIS  ROUTINF  IS  USED  WHEN  INIT  IS  TRUE  TO  GET  STARTED. 

IT  SIMPi.y  COPIES  THE  INITIAL  GUESS  AND  CALLS 

GETCO. 

()  I MENS  I ON  X ( 1 ) ,RHS  ( 1 ) * W ( 259  « S ) » WK  ( 1 ) 

NP  1 =N  ♦ 1 
DO  10  1 = 1,  -gP  1 
1 0 OHS  ( I ♦ 1 ) =w ( I . 1 ) 

CALL  GETCO (W.RHS.WK ,N,H,  0) 

RF TURN 
END 


SUBROUTINE  GFTCOtW.RHS.WK.N.H.M 

THIS  ROUTINF  TAKES  AN  APPROXIMATION  TO  Y AT  N*1  GRID  POINTS. 
GETS  A NUMERICAL  APPROXIMATION  TO  THE  DERIVATIVE 
EXPANSION  EOR  THE  SPLINE  WHICH  THIS  METHOD  SEEKS  AND  THEN 
CALLS  LEQT1B  TO  SOLVE  ThE  LINEAR  EQUATIONS  FOR  THF  SPLINE 
COFFF ICIFmTs. 

DIMENSION  X ( 1) . RHS ( 1 ) .W ( 259,5) .WK ( 1) .C (29) , A (29) ,CTE (?9) 

DATA  C/2^0. .2. ,0. ,-2. ,0, ,2. .0, . 1 1 . 3333333333333.0. .-198, ,0. . 
1 2157.  ,0.. -12131. 33  3 33  3 333  3, 0. .- 3.4 822S6666666bb6E5. 0.  . 
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2 19492202. . 0.  .-S. 94299 l48666bbbF 8,0. . 7. 7626941 3533333E9.0.. 

3 6.07798SS4/31 333E 1 1 , 0 . , -6 . 8422 1 476827 75 7E 1 3 , 0 . , 

4 4.03224739307208F15/ 

NP  1 =N*  1 

N P 3 = N ♦ 3 

I F < *<  .GT.  6 ) 9 TOP  2 7 
KK=4» ( K ♦ ) ) 

00  10  1 r 1 ,KK 

10  CTE ( I ) =C ( I > / (H*H) 

CALL  COEGEN (KK ♦ 1 . . A .CTE ) 

Xl=o. 

x 2 = o . 

no  20  i=i . kk 

X 1 = x 1 ♦ A ( 1 ) *4mS (1*1) 

20  X 2 = X 2 ♦ A ( I )cRlJS(NP3-I  ) 

RHS ( 1 ) = X 1 

RhS  (*,P3)  = X 2 

w ( 1 . 3) = w <1 .5) =1 ./ ( h«h) 

W ( 1 .4) =-2 . / (H*H) 

M (N  ♦ 3, 1 ) =W (N*  3 ♦ 3 ) =1 ./ (H«H) 

* ( N ♦ 3 . 2 ) --2./  ( HaH  ) 

NP2=N*2 

00  100  1=2. NP2 
W(I.1)=W(I»5)=0. 

W ( f «?) =W (I .4) = 1 ,/6. 

100  W ( I . 3) =2. /3. 

CALL  LE  :T 16 ( W,N*3.2.2«RHS.Wk , 0 . IER ) 

IFIIFR  . P.  0 ) 9T0P33 
Rh  TUR'-' 

Em 


SUBROUT  INF  BC (ALFO.BVO.GO.ALFN.BvN.GN.H.N.RHS.RHSTOT) 

This  routine  HANDLES  the  boundary  conditions  for  the  case  when 
yp  is  present  in  the  BOUNDARY  conditions,  it  computes  the 

AMOUNT  RY  WHICH  THE  CURRENT  SOLUTION  FAILS  TO  SATISFY  THE 
BOUNDARY  CONDITIONS  AND  RETURNS  THOSE  VALUES  AS  RHS ( 1 ) AND 

rhs ( n ♦ 3 > . this  routine  works  directly  on  the  spline 

COEFFICIENTS  AND  DOES  NOT  TRY  TO  DO  DIFFERENCE  CORRECTION 
ON  YP  AS  IS  REQUIRED  BY  THE  THEORY.  IT  WOULD  PROBABLY  BE 
EASIER  TO  DO  THAT  DIRECTLY  ON  THE  Y AND  YP  VALUES  THAN  ON  THE 
COEFFICIENTS. 

DIMENSION  RHS ( 1 ) .RHS TOT ( 1 ) ,RCOO (3) »RCON(3) 

BCOO ( 1 ) = ( ALF0»H-3.*RV0> / (6.*H> 

BCOO  (2)  =2.<*A|  FO/3. 

BCOO (3)=(ALF0*H»3.*BV0)/(6.*H) 

BCON (1)=(ALFN*H-3.*RVN)/(6.*H) 

BCON (2) =2. »Al  FN/3. 

BCON ( 3) = ( ALFN*H*3.*RVN) / <6.*H) 

Tfc  1 =0  . 

TE?=0. 

DO  10  1=1.3 

TF 1 =TF 1 .BCOO ( I ) “RHSTOT ( I ) 


L 


n r>  n 


DATA  N‘'A<m  /\>H/ 

CHECK  T HR  VALIOITV  OE  THE  INPUT  PARAMETERS 


SPLIDC 

10  TE2  = IF2.WC0N  ( I ) *RhSTOT (N» I ) 

»HS  ( 1 ) = G 9 - T E 1 
RHS ( N ♦ 1 ) =GN-T)  ? 

RP  TURN 

END 


SUBROUT  INE  L INST (GO. ALFO . UN. ALFN. A .8«N.H. RHSTOT) 

0 I MENS  I ON  RhSTOT  i 1 ) 

THIS  ROUTINE  PRODUCES  The  SPi  INE  COE ~F I C I ENTS  FOR  A 
LINEAR  INTERPOLATION  OE  THE  BOUNDARY  DATA  WHEN  YP  IS  NOT 
INVOLVED  IN  ThE  BOUNDARY  CONDITIONS. 

Y 1 =G0  / A|  rn 
YN  = GN/Ai  N 
SLOP= ( YN-Y  1 ) / ( B- A ) 

RhSTOT (2) =Y1 
DO  30  I - 1 , N 

30  RHSTOT ( I .?) rRHSTOT ( 2 ) ♦ I *H*Sl OP 
RHSTOT ( 1 ) =2.*PHSTOT (2) -RhSTOT l 3) 

RhSTOT (n*3)=2.*RHST0T (N*2) -RHSTOT  < N ♦ 1 ) 

RE' TURN 

END 


c 

c 

C 

C 

C 

C 

C 


C- 

LEQT IB 

— c — — — 

--L  1 

[RRARY  3 

c 

function 

_ 

MATRIX  DECOMPOSITION.  LINEAR  EQUATION 

c 

SOLUTION  - SPACE  ECONOMISE*  SOLUTION  - 

c 

HAND  STORAGE  MODE 

c 

USAGE 

- 

CALL  LEOT1B  ( A , N . NL C . NUC . B . I JOB , I FR ) 

c 

PARAME  TERS 

A 

- 

INPUT/OUTPUT  MATRIX  OF  DIMENSION  N BY 

c 

(NUC*NLC*1).  SEE  PARAMETER  IJOB. 

c 

N 

- 

ORDFR  OF  MATRIX  A AND  THE  NUMBER  OF  ROWS  IN 

c 

H.  (INPUT) 

c 

NLC 

- 

NUMBER  OE  LOWER  CODIAGONALS  IN  MATRIX  A. 

c 

( INPUT ) 

c 

NUC 

- 

NUMBER  OE  UPPER  COOIAGONALS  IN  MATRIX  A. 

c 

( INPUT ) 

c 

I A 

- 

WOW  DIMENSION  OE  A AS  SPECIrlFO  IN  THE 

c 

CALLING  PROGPAM.  (259) 

c 

H 

- 

INPUT/OUTPUT  MATRIX  OE  DIMENSION  N BY  M. 

SUBROUTINE  LEOT1B  ( A . N . NLC . NUC . B , XL . 1 JOB ♦ I ER ) 

THIS  ROUTINE  IS  ADAPTED  FROM  THE  ROUTINE  OF  THE  SAME  NAME 
IN  THE  I MSL  LIBRARY.  THE  CHANGES  ARE: 

(1)  THE  ROW  DIMENSION  OF  A AND  XL  IS  FIXED  AT  259. 

(?)  R IS  A VECTOR  RATHER  THAN  A MATRIX. 


ooo  oooooooonnooooooooonnooonononoonoo 


C 


c 


- • * • M w I y I ^ 

FILL  IN)  Y VALUES  AF  Tf  R DOUBLING  N 
CALL  r,piVAL  (N.H.RHST0T.Y.YP.AY2P) 

WE  USE  DIFFERENT  YP  VALUES  AFTER  AN 


INTERPOL  AT  ION 


C 

C 


S 


XL 

IH 

I JOB 


IFR 


DIMENSION 

DATA 

IFR  = 0 

JBEG  = NL  C ♦ 1 

NLC 1 = (“EG 

IF  fIJOR  .EG.  2) 

RN  = N 


! = 1 

NC  = JBEG*NUC 
NN  = NT 
JEND  = N' 

IF  (N  .EG.  1 .OR 
k = 1 


ON  INPUT.  R CONTAINS  THE  M RlOHT-HAND  SIDES 
OF  THE  EQUATION  AX  = B.  ON  OUTPUT,  THE 
SOLUTION  MATRIX  X REPLACES  B.  IE  IJOB  = 1, 

H IS  NOT  USED. 

THIS  VERSION  ASSUMES  ONE  RIGHT  HAND  SIDE 
INPUT/OUTPUT  MATRIX  OF  SIZE  777  WHICH  CONTAINS 
PIVOTING  INFORMATION  WHICH  SHOULD  NOT  BE 
DESTROYED  BY  THE  CALLING  PROGRAM. 

ROW  DIMENSION  OF  B ^S  SPECIFIED  IN  THE 
CALLING  PROGRAM.  (259) 


- INPUT  OPTION  PARAMETER.  I JOB  = I IMPLIES  WHEN 
I = 0,  FACTOR  THE  MATRIX  A AND  -OLVE  THE 
EQUATION  AX  = B.  ON  INPUT,  A CONTAINS  THE 
COEFFICIENT  MATRIX  OF  T HE  EQUATION  AX  = B. 
WHERE  A IS  ASSUMED  TO  BE  AN  N BY  N BAND 
MATRIX.  A IS  STORED  IN  RAND  STORAGE  MODE 
AND  THEREFORE  HAS  DIMENSION  N BY 
( NLC *NUC ♦ 1 ) • ON  OUTPJT , A IS  REPLACEO 
BY  THE  U MATRIX  OF  THE  L-U  OEC0MP0S I T I ON 
OF  A ROWWISE  PERMUTATION  OF  MATRIX  A.  U IS 
STORED  IN  BAND  STORAGE  mODE. 

1=1,  FACTOR  THE  matrix  A.  A CONTAINS  THE 


SAME  INPUT/OUTPUT  INFORMATION  AS  IE 
IJOB  = 0. 

I = 2.  SOI  VF  THE  EQUATION  AX  = B.  THIS 

option  implies  that  leotib  has  already 

BEEN  CALLFD  USING  IJOB  = 0 OR  1 SO  THAT 
THE  MATRIX  A HA'.  ALREADY  BEEN  FACTORED. 
IN  THIS  CASE.  OUTPUT  MATRIX  A 
'•'US T have  BEEN  SAVED  FOR  REUSE  IN  THE 
CALL  TO  LFQT1B. 

- ERROR  PARAMETER. 

TERMINAL  ERROR  = 

N = 1 INDICATES  THAT  MATRIX  A IS 

ALGORITHMICALLY  SINGULAR.  (SEE  THE 
CHAPTER  L PRELUDE) . 

A (259,1) .XL (259,3) ,B(259> 

7ERO/0./, ONE/1 .0/ 


GO  TO  «0 

RESTRUCTURE  THE  ■ A TR I X 

FIND  RECIPROCAL  OF  THE  LARGE  9 T 

4BS0FUTE  VALUE  IN  ROW  I 


NLC  .EQ.  0)  GO  TO  25 


P = ZE Ro 

DO  10  J = JBEG.JEND 
A ( I ♦ K ) = A ( I « J ) 

Q = 4RS ( A ( t ,K  ) | 

IF  (0  .GT.  P)  P = 0 


I s } I 1-tt  MU“Ht  K U*  t I I L 1 " SU"iVt  "UCT'U.c 

IS  NOT  EXPECTED  MV  ITLIMM. 

IE  THE  Change  IN  v (RNORmi  INCREASES  BUT  THE  RESIDUAL  DECREASES 


74 


SPLIDC  - ?5 


K = K » 1 

10  CONTINUE 

IE  IP  .FO.  ZERO)  GO  TO  135 
XL ( I » NL  C 1 ) = ONE/P 
IE  (K  .gt.  NO  GO  TO  20 
DO  15  ) = K.NC 
A ( I * ) ) = ZERO 

15  CONTINUE 
20  I = 1*1 

JHEG  = JHEG-1 

IE  (JEND-JBEG  .EO.  N)  JLND  = JEND-1 
IE  (I  .LF . NEC ) GO  TO  5 
JHEG  = I 
NN  = JEN' 

?5  JENO  = N-NUC 

DO  SO  I JBE G . N 
P = TCRO 
DO  30  ) = 1 . NN 

0 = APS ( A ( I , J)  ) 

IE  (0  .GT.  P)  P r Q 
30  CONTINUE 

IE  (P  .EO.  ZERO)  GO  to  135 
XL  l I »NLC1  ) = ONE/P 
IE  (I  .EO.  JEND)  GO  TO  37 
IE  (I  ,LT.  JEND)  GO  TO  SO 
x - NN ♦ 1 
DO  )5  J = K.NC 
A(I,J)  = /PRO 
35  CONTINUE 

37  NN  = NN-1 

SO  CONTINUE 
L = NL  C 

L-U  DECOMPOSITION 

DO  75  k-  - i,n 

P = ABS ( A (K* 1 ) >*XL (K.NLC1 ) 

T = x 

IE  (L  .LT.  N)  L = L*1 
k 1 = k ♦ 1 

IP  (-1  .GT . L ) GO  TO  50 
DO  S5  i = K 1 . L 

= ABS ( A ( J. 1 ) ) »XL ( J.NLC1 ) 

IE  (Q  ,LE.  P)  GO  TO  S5 
p 0 
T - J 
s5  CONTINUE 

50  XL ( I.NLC1 ) = XL (K.NLC1  ) 

XL(K.NLCI)  = I 

slNGUI  ARITY  FOUND 
IE  <Rv*P  .EO.  RN)  GO  TO  135 

INTERCHANGE  ROWS  I AND  < 

IE  («■  .EO.  I)  GO  TO  SO 

DO  55  i = 1 ,NC 
P s A (K  , J) 

A ( K « J ) = A ( I * J ) 

A (I . J)  = P 
55  CONTINUE 

SO  IE  < k 1 .GT.  L)  GO  TO  75 


o o n 


POSSIBLE.  IF  THE  FINAL  nfwton  iteratf  solves  the  non- 
l INFAR  SYSTEM  better  than  the  FIRST  ITERATE,  then  we  get 
an  FRROR  ESTIMATE  FOR  USE  IN  TESTING  ORDER  AND  WE  USE 


SPLIDC  -2b  < 

o 70  T = X 1 .L 

P - A ( I . 1 ) / A ('  .1) 

IK  = I-K 

<L  (K  I . IK  ) = P 

DO  b5  J = 2,  NO 

A ( I . J-  1 ) = A ( I . J) -P»A  <K , J ) 

b5  CONTINUE 

A(I.NC)  = 7ER0 
70  CONTINUE 

75  CONTINUE 

IF  (I  JOB  ,E<' . 1 ) GO  TO  90  05 

FORW'RD  SUBSTITUTION 

SOL  = NLC 


105 

H 

= 

1 

.N 

I = 

YL 

( 

* ♦ 

NLC1  ) 

IF 

( I 

• 

F 0 

. K) 

GO 

TO  BO 

p = 

B (K  ) 

B { x 

) 

= 

R ( I ) 

B ( I 

) 

= 

P 

IF 

(L 

• 

LT 

. N) 

L = 

L ♦ 1 

K 1 

= K 

♦ 

1 

IF 

(K  1 

.GT.  L) 

GO 

TO  105 

DO 

100 

I 

= * 1 * 

L 

IK 

= 

I 

-K 

P = x L ( k 1 , I K ) 


B ( I ) = B(I) -P*B ( K ) 

100  CONTINUE 
105  CONTINUE 

C BACKWARD  SUBSTITUTION 

JBEG  = N C ♦ NL C 
i = 1 
X 1 - N ♦ 1 
rO  120  I - l.N 
K - K 1 -I 
P = B ( K ) 

IF  (L  .FQ.  1)  GO  TO  115 
no  110  KK  = 2 «L 


IK  = KK  ♦ K 

P = P-A (K .KK ) »R ( IK-  1 ) 


110 

CONTINUE 

115 

B (K)  = 

P/A (K , 1 ) 

IF  (L 

.LE.  JBEG)  L = L ♦ 1 

120 

CONTINUE 

GO  TO  9005 

135 

IF  P = 12  ) 

9000 

CONI  INUF 

9005 

return 

END 


POIMT  a 70 » I NOE  * tPFSNtPNOW^ 

CALL  PWSUB (N*K 1 .ERRNIOP. X .NhSTOT ) 
CONTINUE 


100  1 
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The  methodology  is  very  briefly  described  and  then  numerical  results  are 
presented  for  an  implementation  of  a smooth-cubic-spline-collocation  procedure 
accelerated  via  iterated  deferred  corrections  to  obtain  approximate  solutions, 
accurate  to  a prescribed  tolerance,  of  two-point  boundary-value  problems  for 
second-order  scalar  ordinary  differential  equations.  The  results  are  similar 
to  those  obtained  via  a less  generally  applicable  finite-difference-oriented 
code  described  elsewhere  which  competes  well  with  the  best  codes  available. 
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